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ABSTRACT 


Airspace  encounter  models,  covering  close  encounter  situations  that  may  occur 
after  standard  separation  assurance  has  been  lost,  are  a  critical  component  in  the 
safety  assessment  of  aviation  procedures  and  collision  avoidance  systems.  Of  par¬ 
ticular  relevance  to  Unmanned  Aircraft  Systems  is  the  potential  for  encountering 
general  aviation  aircraft  that  are  flying  under  Visual  Flight  Rules  (VFR)  and  which 
may  not  be  in  contact  with  air  traffic  control.  In  response  to  the  need  to  develop 
a  model  of  these  types  of  encounters,  Lincoln  Laboratory  undertook  an  extensive 
radar  data  collection  and  modeling  effort.  This  report  describes  the  structure  and 
content  of  that  encounter  model.  Bayesian  networks  are  used  to  represent  relation¬ 
ships  between  dynamic  variables  and  to  construct  random  aircraft  trajectories  that 
are  statistically  similar  to  those  observed  in  the  radar  data. 

The  model  described  in  this  report  is  one  of  three  that  are  currently  under 
development  by  Lincoln  Laboratory.  An  encounter  with  an  intruder  that,  does  not 
have  a  transponder,  or  between  two  aircraft  using  a  Mode  A  code  of  1200  (VFR), 
is  uncorrelated  in  the  sense  that,  it  is  unlikely  that  there  would  be  prior  interven¬ 
tion  by  air  traffic  control.  The  uncorrelated  model  described  in  this  report  is  based 
on  transponder-equipped  aircraft  using  a  1200  (VFR)  Mode  A  code  observed  by 
radars  across  the  U.S.  In  addition  to  representing  VFR-on-VFR  encounters,  this 
model  is  representative  of  encounters  between  a  cooperative  aircraft  and  conven¬ 
tional  non-cooperative  aircraft  similar  to  those  that  use  a  VFR  transponder  code. 
A  second  uncorrelated  model  is  also  being  developed  based  on  true  primary-only 
(non-cooperative)  aircraft  tracks,  though  the  quality  of  such  tracks  is  inferior  to 
that  available  from  the  cooperative  1200-code  tracks  used  here.  Finally,  a  corre¬ 
lated  encounter  model  is  also  under  development  to  represent  situations  in  which  it 
is  likely  that  there  would  be  air  traffic  control  intervention  prior  to  a  close  encounter. 
The  correlated  model  will  apply  to  aircraft  that  are  using  a  discrete  (non- 1200)  mode 
A  code. 

A  separate  electronic  file  is  provided  by  Lincoln  Laboratory  that  contains 
the  statistical  data  required  to  generate  encounter  trajectories.  Details  on  how  to 
interpret  the  data  file  and  an  example  of  randomly  constructing  trajectories  are 
provided  in  Appendices  A  and  B,  respectively. 
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1.  INTRODUCTION 


One  of  the  main  challenges  to  integrating  unmanned  aircraft  into  the  National  Airspace  Sys¬ 
tem  is  the  development  of  systems  that  are  able  to  sense  and  avoid  local  air  traffic.  If  designed 
properly,  these  collision  avoidance  systems  could  provide  an  additional  layer  of  protection  that 
maintains  or  even  enhances  the  current  exceptional  level  of  aviation  safety.  However,  due  to  their 
safety-critical  nature,  rigorous  assessment  is  required  before  sufficient  confidence  can  exist  to  certify 
collision  avoidance  systems  for  operational  use.  Evaluations  typically  include  flight  tests,  opera¬ 
tional  impact  studies,  and  simulation  of  millions  of  traffic  encounters  with  the  goal  of  exploring 
the  robustness  of  the  collision  avoidance  system.  Key  to  these  simulations  are  so-called  encounter 
models  that  describe  the  statistical  makeup  of  the  encounters  in  a  way  that  represents  what  actually 
occurs  in  the  airspace. 

One  system  that  has  been  rigorously  tested  in  this  manner  is  the  Traffic  alert  and  Collision 
Avoidance  System  (TCAS).  As  part  of  the  TCAS  certification  process  in  the  1980s  and  1990s, 
several  organizations  tested  the  system  across  millions  of  simulated  close  encounters  and  evaluated 
the  risk  of  a  near  mid-air  collision  (NMAC,  defined  as  separation  less  than  500  ft  horizontally 
and  100  ft  vertically)  [1-4].  This  analysis  ultimately  led  to  the  certification  and  U.S.  mandate 
for  TCAS  equipage  on  larger  transport  aircraft.  More  recently,  Eurocontrol  and  the  International 
Civil  Aviation  Organization  (ICAO)  performed  similar  sets  of  simulation  studies  for  European  and 
worldwide  TCAS  mandates  [5,6]. 

The  design  of  a  collision  avoidance  system  represents  a  careful  balance  to  enhance  safety 
while  ensuring  a  low  rate  of  unnecessary  maneuvers.  This  balance  is  strongly  affected  by  the  types 
of  encounter  situations  to  which  the  system  is  exposed.  It  is  therefore  important  that  simulated 
encounters  are  representative  of  those  that  occur  in  the  airspace.  Hence,  tremendous  effort  has 
been  made  by  various  institutions  since  the  early  1980s  to  develop  encounter  models  based  on 
radar  data  [1,3,7  10].  The  primary  contribution  of  this  report  is  to  introduce  a  new  approach  to 
encounter  modeling  that  is  based  on  a  Bayesian  statistical  framework  and  which  uses  recent  radar 
data  collected  across  the  U.S.  The  advantage  of  applying  a  Bayesian  framework  is  that  it  allows 
us  to  optimally  leverage  available  radar  data  to  produce  a  model  that  is  representative  of  actual 
operations. 


1.1  ENCOUNTER  TYPES 

The  encounters  covered  by  this  model  are  those  involving  aircraft  in  the  final  stages  before 
a  collision  over  a  period  of  time  on  the  order  of  one  minute.  It  is  assumed  that  prior  safety  layers 
(e.g.,  airspace  structure,  Air  Traffic  Control  (ATC)  advisories  or  vectors)  have  failed  to  maintain 
standard  separation  distances  between  aircraft.  The  model  is  therefore  useful  in  describing  the 
types  of  situations  that  need  to  be  addressed  by  a  collision  avoidance  system,  but  will  not  address 
earlier  aspects  such  as  ATC  communications  or  coordination. 

There  are  two  fundamental  types  of  close  encounters.  In  the  first,  both  aircraft  involved  are 
cooperative  (i.e.,  have  a  transponder)  and  at  least  one  is  in  contact  with  air  traffic  control.  It  is  then 
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likely  that  at  least  one  aircraft  will  receive  some  notification  about  the  traffic  conflict  and  begin 
to  take  action  before  a  collision  avoidance  system  gets  involved.  We  term  this  type  of  encounter 
“correlated”  because  the  trajectories  of  each  aircraft  may  involve  maneuvers  that  are  correlated  to 
some  degree  due  to  this  prior  intervention.  The  second  type  of  encounter  we  term  “uncorrelated” 
and  involves  at  least  one  noncooperative  aircraft  or  two  aircraft  flying  under  Visual  Flight  Rules 
(VFR)  without  flight  following  (i.e.,  using  a  transponder  Mode  A  code  of  1200).  In  these  encounters, 
it  is  unlikely  that  air  traffic  control  would  become  involved  prior  to  the  close  encounter;  rather  the 
two  aircraft  must  rely  solely  on  visual  acquisition  (or  some  other  collision  avoidance  system)  at  close 
range  to  remain  separated.  Such  encounters  tend  to  be  uncorrelated  since  there  is  no  coordinated 
intervention  prior  to  the  close  encounter.  The  assumption  behind  uncorrelated  encounters  is  that 
the  two  aircraft  blunder  into  close  proximity;  it  is  up  to  the  collision  avoidance  system  not  the 
encounter  model  -to  describe  how  the  aircraft  respond  to  the  encounter.  A  complete  evaluation  of 
unmanned  systems  will  require  analysis  using  both  correlated  and  uncorrelated  models. 

This  report  focuses  on  uncorrelated  encounter  modeling.  In  this  report,  we  base  our  uncorre¬ 
lated  model  on  beacon  reports  from  aircraft  squawking  1200,  not  radar  returns  from  noncooperative 
traffic.  Radar  surveillance  of  noncooperative  targets  is  typically  complicated  due  to  clutter  and 
missed  detections,  making  identification  of  real  tracks  difficult.  The  lack  of  a  transponder  means 
that  only  position  information,  and  no  identity  code  or  altitude,  is  available.1  Hence,  it  is  difficult 
to  infer  vertical  rates,  an  important  component  of  an  encounter  model. 

Beacon-equipped  aircraft  can  transmit  either  a  discrete  Mode  A  code  or  the  non-discrete 
code  of  1200.  Aircraft  using  code  1200  tend  to  be  general  aviation  typically  flying  VFR.  They 
generally  fly  at  low  altitudes  and  make  significantly  more  maneuvers  than  transport  aircraft,  both 
horizontally  and  vertically.  Thus  to  a  large  degree  they  resemble  noncooperative  aircraft.  Due 
to  the  poor  quality  of  noncooperative  data,  we  use  1200  tracks  here  as  surrogates  for  primary- 
only  tracks  for  the  construction  of  this  model.  Other  work  is  underway  toward  processing  true 
non-cooperative  tracks  and  will  be  described  in  a  future  document. 

The  1200  tracks  are  a  suitable  surrogate  for  much  of  the  noncooperative  traffic  in  the  National 
Airspace  System,  but  there  are  certain  categories  of  noncooperative  targets  for  which  they  are  not 
suitable.  For  example,  most  balloons,  ultralights,  and  gliders  do  not  fly  like  transponder-equipped 
aircraft  squawking  1200.  The  challenge  in  developing  models  for  unconventional  aircraft  such  as 
balloons,  ultralights,  and  gliders  is  the  lack  of  high-quality  data.  As  these  data  become  available, 
the  same  processes  outlined  in  this  report  will  be  used  to  estimate  a  noncooperative  model  for 
unconventional  aircraft.  Additionally,  where  such  data  are  not  available,  a  model  can  be  manually 
constructed  based  on  knowledge  of  typical  flight  trajectories  and  performance  limits.  Table  1  shows 
which  encounter  model  to  use  depending  on  the  types  of  aircraft  involved  in  the  encounter. 


^ome  military  radars  have  height-finding  capability,  although  the  accuracy  of  the  altitude  generated  is  significantly 
inferior  to  the  transponder  Mode  C  altitude. 
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TABLE  1.  Encounter  model  types.  There  are  three  types  of  encounter  models:  correlated  (C),  uncorrelated 
conventional  VFR  (V),  and  uncorrelated  unconventional  VFR  (U).  This  table  shows  which  model  to  use  depending 
on  the  types  of  aircraft  involved  in  the  encounter.  In  terms  of  aircraft  types,  Discrete  indicates  a  cooperative 
aircraft  using  a  non-1200  Mode  A  transponder  code;  VFR  denotes  a  cooperative  aircraft  using  the  1200  Mode  A 
transponder  code;  Non-coop.  Conv.  indicates  a  non-cooperative  but  conventional  fixed-wing  aircraft;  Non-coop 
Unconv.  indicates  a  non-cooperative  and  unconventional  aircraft  (e  g.,  balloon,  ultralight).  Only  the  uncorrelated 
conventional  VFR  model  (in  bold)  is  discussed  in  this  report;  the  other  models  are  under  development. 


Intruder  Aircraft 

Aircraft  of  Interest 

Discrete 

VFR 

Discrete 

c 

C 

VFR 

c 

V 

Non-coop.  Conv. 

V 

V 

Non-coop.  Unconv. 

u 

u 

1.2  PROCESS  OVERVIEW 

Figure  1  diagrams  the  steps  involved  in  processing  radar  data  to  build  the  encounter  model. 
The  high-level  approach  is  to  model  nominal  aircraft  trajectories  using  Markov  models  represented 
by  Bayesian  networks  (Section  2).  The  first  step  in  constructing  a  Markov  model  involves  extract¬ 
ing  features,  such  as  turn  rate  and  vertical  rate,  from  the  radar  data.  From  these  features,  we 
collect  sufficient  statistics  that  describe  the  distribution  over  maneuvers  and  other  properties  of 
trajectories.  We  use  Bayesian  model  selection  methods  to  search  for  the  best  network  structure  for 
representing  the  observed  data.  (Section  3).  We  use  the  best  network  structure  and  the  associated 
sufficient  statistics  to  generate  new  trajectories  that  are  representative  of  those  observed  by  radar 
(Section  4).  The  trajectories  are  then  used  to  simulate  encounters  between  aircraft  (Section  5).  A 
Monte  Carlo  safety  analysis  involves  running  a  large  collection  of  encounters  and  collecting  data 
about  near  mid-air  collision  (NMAC)  frequency  (Section  6). The  final  section  of  this  report  summa¬ 
rizes  and  discusses  further  work. 

The  main  outputs  of  the  work  described  in  this  report  are  the  sufficient  statistics  and  model 
structure.  A  digital  representation  of  the  sufficient  statistics  and  model  structure  are  publicly 
available  from  Lincoln  Laboratory  (Appendix  A).  An  example  of  using  the  data  in  that  file  to 
construct  a  random  trajectory  is  provided  in  Appendix  B. 

Throughout  this  report  we  use  the  standard  units  used  in  aviation.  In  particular,  altitudes 
are  in  feet,  positions  are  in  nautical  mile  coordinates,  vertical  rates  are  in  feet  per  minute,  turn 
rates  are  in  degrees  per  second,  airspeeds  are  in  knots  (true  airspeed),  and  accelerations  are  in 
knots  per  second.  Time  is  reported  in  seconds. 
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Modeling 


Simulation 


FIGURE  1.  Modeling  and  simulation  process  overview.  The  sufficient  statistics  and  network  structure  (in  bold) 
are  the  main  elements  provided  as  part  of  this  work. 


4 


2.  MODEL 


We  model  nominal  flight,  i.e.  flight  without  avoidance  maneuvering,  using  a  Markov  process 
represented  by  a  dynamic  Bayesian  network.  A  Markov  process  is  a  stochastic  process  where  the 
probability  distribution  over  future  states  is  conditionally  independent  of  past  states  given  the 
present  state.  In  other  words,  one  only  needs  to  know  the  present  state  to  predict  the  next  state. 

The  states  in  our  model  specify  how  the  position,  altitude,  and  airspeed  change  over  time. 
In  particular,  each  state  specifies  a  vertical  rate  /*,  turn  rate  i/>,  and  airspeed  acceleration  v.  Given 
an  initial  airspeed  v,  horizontal  coordinates  ( x,y ),  heading  xj),  vertical  rate  h ,  altitude  layer  L,  and 
airspace  class  A,  we  can  infer  from  our  model  how  the  aircraft  trajectory  evolves  over  time. 

One  way  to  represent  a  Markov  model  is  with  an  exhaustive  state-transition  matrix  that 
specifies  the  probability  of  transitioning  between  all  pairs  of  states.  Unfortunately,  the  number  of 
independent  parameters  required  to  define  the  matrix  grows  super  exponentially  with  the  number 
of  variables  defining  the  model.  The  more  independent  parameters  there  are  in  the  model,  the  more 
data  one  needs  to  properly  estimate  their  values.  However,  by  using  dynamic  Bayesian  networks, 
we  can  leverage  conditional  independence  between  some  variables  to  greatly  reduce  the  number 
of  parameters.  We  can  learn  the  structure  of  the  dynamic  Bayesian  network  by  maximizing  the 
posterior  probability  of  the  network  structure  given  the  data.  Appendix  A  provides  the  necessary 
background  on  Bayesian  networks.  The  remainder  of  this  section  outlines  our  approach. 


2.1  MODEL  VARIABLES 

There  are  six  variables  in  the  uncorrelated  encounter  model: 

•  Airspace  class  A:  This  variable  may  take  on  one  of  four  values:  B,  C,  D,  and  O,  indicating 
which  class  of  airspace  the  aircraft  is  in.  The  values  B,  C,  and  D  correspond  to  the  controlled 
airspace  classes  defined  by  the  FAA.  The  value  O  represents  “other  airspace,”  that  is  airspace, 
such  as  Class  A,  E,  G,  that  is  not  B,  C,  or  D.  The  airspace  class  variable  was  incorporated 
into  our  model  to  account  for  the  variation  in  how  aircraft  fly  in  different  airspace  classes. 
Note  that  there  should  be  no  VFR  aircraft  in  Class  A  due  to  the  requirement  that  aircraft  in 
that  Class  of  airspace  fly  under  Instrument  Flight  Rules. 

•  Altitude  layer  L :  Airspace  is  also  divided  into  four  altitude  layers,  in  a  process  similar  to 
prior  encounter  models  developed  by  Eurocontrol.  The  first  layer  spans  from  500  to  1200  ft 
Above  Ground  Level  (AGL)  to  capture  aircraft  in  the  traffic  pattern  or  performing  low- 
level  maneuvers.  The  second  layer  spans  a  transition  zone  from  1200  to  3000  ft  AGL,  the 
cruise  altitude  where  the  hemispheric  rule  begins.  The  third  layer  spans  from  3000  ft  AGL 
to  5000  ft  AGL  covering  a  mix  of  low-altitude  enroute  and  maneuvering  aircraft.  The  fourth 
layer  includes  airspace  above  5000  ft  AGL  and  would  cover  most  enroute  VFR  traffic. 

•  Airspeed  v:  We  model  true  airspeed  and  allow  it  to  vary  during  flight. 
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•  Acceleration  v:  Unlike  previous  encounter  models,  we  allow  airspeed  acceleration  to  vary 
at  every  second. 

•  Turn  rate  ip:  Turn  rate  is  permitted  to  change  every  second  in  our  model.  The  prior 
European  and  ICAO  cooperative  models  allowed  only  a  single  turn  during  an  encounter. 

•  Vertical  rate  h:  The  vertical  rate  is  permitted  to  change  at  every  second.  All  prior  cooper¬ 
ative  models  allowed  only  a  single  vertical  acceleration  period  during  an  encounter. 

The  remainder  of  this  section  explains  how  we  model  the  joint  probability  distribution  over  these 
variables.  A  core  concept  of  the  model  is  that  first,  parameters  are  randomly  selected  to  describe 
the  situation  at  the  beginning  of  an  encounter  simulation.  Second,  certain  parameters  may  vary 
dynamically  during  an  encounter  so  that  realistic  turning,  climbing,  or  accelerating  trajectories 
can  be  generated.  Thus  there  are  two  separate  probability  distributions  in  the  model:  an  initial 
distribution  to  set  up  an  encounter  situation,  and  a  transition  distribution  to  describe  how  the 
trajectory  evolves  over  time. 

The  probability  distributions  are  represented  using  a  Bayesian  network.  The  network  de¬ 
scribes  which  parameters  have  dependencies  on  other  parameters.  For  example,  it  would  be  ex¬ 
pected  that  the  turn  rate  of  an  aircraft  would  be  dependent  on  its  vertical  rate.  Without  properly 
capturing  this  dependency,  unrealistic  situations  may  be  generated,  e.g.,  involving  aircraft  with 
simultaneously  high  climb  rates  and  high  turn  rates.  A  Bayesian  scoring  process  (Appendix  A) 
was  used  to  evaluate  different  Bayesian  network  structures  and  arrive  at  a  structure  that  optimally 
represents  the  observed  radar  data. 


2.2  INITIAL  DISTRIBUTION 

In  our  model  of  aircraft  flight,  we  are  concerned  with  modeling  the  distribution  over  the  initial 
values  of  h,  ip,  v ,  v,  L,  and  A.  Figure  2a  shows  the  structure  we  use  for  our  modeling  effort.  Other 
structures  are  certainly  plausible,  and  it  is  possible  to  compare  different  structures  using  Bayesian 
scoring  (Appendix  C)  to  determine  which  structure  is  more  likely  given  the  data.  We  chose  the 
network  structure  in  Figure  2a  because  it  had  the  highest  score  among  the  selection  of  candidate 
networks  we  considered  (see  Appendix  G). 

Given  a  structure,  sufficient  statistics  extracted  from  data,  and  a  Bayesian  prior,  we  then 
sample  from  the  Bayesian  network  to  produce  initial  airspace  classes,  altitude  layers,  vertical  rates, 
turn  rates,  airspeeds,  and  accelerations  that  are  representative  of  those  found  in  the  data.  The 
boxes  and  arrows  used  in  the  structural  diagram  show  the  order  in  which  this  sampling  occurs.  For 
example,  based  on  the  structure  in  Figure  2a,  to  determine  the  initial  state  of  the  aircraft  we  first 
randomly  determine  an  airspace  class  A.  Once  the  airspace  class  has  been  determined,  an  altitude 
layer  L  is  selected.  The  probabilities  associated  with  each  altitude  layer  depend  on  which  airspace 
class  was  chosen  earlier.  Once  A  and  L  have  been  selected,  we  then  randomly  select  airspeed  v, 
and  so  on.  An  example  of  working  through  this  process  is  provided  in  Appendix  B.  Alternately, 
we  could  assign  outright  an  airspace  class  and/or  altitude  layer  for  a  particular  study  and  then 
randomly  select  values  for  the  remaining  variables. 
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2.3  TRANSITION  DISTRIBUTION 


We  use  a  separate  Bayesian  network  to  model  how  the  variables  h,  ip,  and  v  evolve  over 
time.  The  structure  we  use  is  shown  in  Figure  2b.  In  this  network,  we  have  a  first  layer  that 
represents  the  state  of  the  system  at  the  present  time  step  and  a  second  layer  that  represents  the 
state  of  the  system  at  the  next  time  step.  There  may  be  dependencies  between  layers  and  within 
the  second  layer.  Such  a  two-layer  temporal  Bayesian  network  is  known  as  a  dynamic  Bayesian 
network  [11,12].  Parameter  and  structure  learning  in  dynamic  Bayesian  networks  is  similar  to 
regular  Bayesian  networks  (Appendix  C).  Again,  we  chose  the  highest-scoring  network  structure 
among  our  candidate  network  structures  (see  Appendix  G).  Given  a  structure,  sufficient  statistics 
extracted  from  data,  and  a  prior,  we  then  sample  from  the  Bayesian  network  to  determine  the  next 
vertical  rate,  turn  rate,  and  acceleration  command  that  are  representative  of  what  was  observed  in 
the  data. 

In  general,  time  steps  in  dynamic  Bayesian  networks  may  be  of  any  duration,  but  for  our 
modeling  effort  we  chose  steps  of  1  s.  Shorter  time  steps  allow  for  more  frequent  variations  in 
airspeed,  vertical  rate,  and  turn  rate,  but  they  require  more  computation  per  unit  of  simulation 
time.  Time  steps  of  1  s  balance  maneuver  complexity  with  computation. 

A  complete  trajectory  is  constructed  by  updating  the  aircraft  state  in  1  s  intervals.  Within 
each  interval,  the  three  derivative  variables  h,  ip,  and  v  are  treated  as  target  values  and  held 
constant.  A  dynamic  model  (which  is  beyond  the  scope  of  this  report)  is  used  to  compute  and 
update  the  aircraft  state  at  each  time  step  based  on  these  piecewise-constant  target  values. 
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(a).  Initial  distribution.  (b>.  Transition  distribution 


FIGURE  2.  Bayesian  networks  representing  the  variable  dependency  structure  for  the  initial  and  transition 
distributions. 
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3.  ESTIMATION 


This  section  describes  the  processing  required  to  transform  radar  tracks  into  sufficient  statis¬ 
tics  that  may  be  used  to  model  uncorrelated  encounters.  Figure  3  outlines  this  process. 


3.1  RADAR  TRACKS 

Our  radar  data  stream  comes  from  the  84th  Radar  Evaluation  Squadron  (RADES)  at  Hill 
AFB,  Utah.  RADES  receives  radar  data  from  FAA  and  Department  of  Defense  sites  throughout 
the  United  States.  They  maintain  continuous  real-time  feeds  from  a  network  of  sensors,  including 
long-range  ARSR-4  radars  around  the  perimeter  of  the  United  States  and  short-range  ASR-8,  ASR- 
9,  and  ASR-11  radars  in  the  interior.  Radar  ranges  vary  from  60  to  250  NM.  Figure  4  shows  the 
coverage  by  the  more  than  120  sensors  whose  data  was  used  to  construct  our  model.  Figure  5  shows 
the  density  of  VFR  traffic,  and  Appendix  D  plots  flight  hour  densities  at  different  altitude  layers. 
Appendix  E  summarizes  the  volume  of  data  received  from  each  sensor. 

Note  that  there  are  a  number  of  advantages  to  our  RADES  data  feed  compared  to  the  En¬ 
hanced  Traffic  Management  System  (ETMS)  data  often  used  in  airspace  analyses.  ETMS  data 
include  only  cooperative  aircraft  on  filed  Instrument  Flight  Rules  flight  plans  and  provides  updates 
once  per  minute  showing  aircraft  position  after  processing  by  air  traffic  control  automation.  In  con¬ 
trast,  RADES  data  is  continuously  streaming  directly  from  the  radar,  includes  primary-only  radar 
returns  as  well  as  all  cooperative  transponder  returns  (whether  on  a  flight  plan  or  not),  providing 
track  updates  every  5  or  12  seconds  without  being  affected  by  automation  systems.  This  ensures 
that  our  filters  and  trackers  have  the  best  raw  data  with  which  to  begin  processing. 

National  Offload  Program  (NOP)  data  is  another  source  that  we  could  have  used  to  construct 
our  model.  An  advantage  of  NOP  data  is  the  inclusion  of  flight-plan  and  aircraft-type  information. 
However,  NOP  data  is  post  automation,  like  ETMS.  NOP  does  not  include  data  from  Department 
of  Defense  sensors  and  does  not  have  as  comprehensive  coverage  as  our  RADES  feed. 

To  build  our  model,  we  collected  VFR  (1200-code)  beacon  reports  between  December  1,  2007 
and  December  7,  2007,  amounting  to  30,000,000  reports  representing  78,000  flight  hours.  The  raw 
radar  data  is  first  processed  using  a  tracking  algorithm  developed  at  Lincoln  Laboratory  [13].  A 
fusion  algorithm,  also  developed  at  Lincoln  Laboratory  [14],  then  fuses  tracks  from  multiple  sensors 
to  give  one  global  view  of  all  the  tracks  in  U.S.  airspace  (see  Appendix  F).  We  eliminated  tracks 
that  had  fewer  than  ten  scans.  We  found  that  approximately  ten  scans  are  required  to  accurately 
estimate  the  various  maneuver  rates.  We  also  eliminated  tracks  if  any  of  their  associated  reports 
were  inside  Special  Use  Airspace  whose  boundaries  are  defined  in  the  Digital  Aeronautical  Flight 
Information  File  (DAFIF),  8th  Edition,  managed  by  the  National  Geospatial-Intelligence  Agency 
(NGA). 
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FIGURE  3.  Estimation  process  flow. 
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FIGURE  4.  Radar  coverage  map. 


FIGURE  5.  VFR  density  between  500ft  AGL  and  FL180  in  log10  aircraft  per  NM3.  The  locations  of  the  sensors 
are  indicated  by  magenta  circles. 
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3.2  OUTLIER  REMOVAL 


The  first  step  in  processing  the  raw  radar  tracks  is  to  detect  and  remove  outliers.  In  the 
horizontal  plane,  we  remove  jumps  with  ground  speeds  above  600  kt  using  the  following  algorithm. 
We  begin  by  estimating  the  speed  between  each  sample  point  by  dividing  the  distance  between 
samples  by  the  time  interval  between  samples.  Samples  on  either  side  of  segments  where  the  speed 
is  above  the  threshold  of  600  kt  are  stored  in  a  list  of  candidates  for  removal.  We  iterate  through 
the  list  of  candidates  and  remove  the  one  that  minimizes  the  sum  of  speeds  above  the  set  threshold. 
The  process  repeats  until  there  are  no  longer  any  segments  with  speeds  above  the  threshold. 

In  the  vertical  plane,  we  remove  missing  Mode  C  altitude  reports  from  consideration.  Then, 
using  the  same  process  as  used  for  the  horizontal  plane,  we  remove  outliers  with  vertical  rates 
greater  than  5000ft/min  or  less  than  —5000  ft /min.  We  remove  altitude  reports  that  come  before 
the  first  position  report  or  after  the  last  position  report  to  prevent  extrapolation.  We  also  remove 
position  reports  that  come  before  the  first  altitude  report  or  after  the  last  altitude  report,  also  to 
prevent  extrapolation.  After  outlier  removal,  we  discard  tracks  with  fewer  than  ten  valid  scans. 


3.3  TRACK  SMOOTHING 

After  removing  any  outliers  from  a  track,  we  smooth  the  remaining  data  points,  first  hori¬ 
zontally  and  then  vertically.  We  use  the  same  smoothing  scheme  for  both  horizontal  and  vertical 
smoothing.  We  use  the  following  general  formula  to  transform  a  raw  trajectory  (t\,  X\), . . . ,  (tn,  xn) 
to  a  smoothed  trajectory  yx, . . . ,  yn. 


T,jW(ti,tj)Xi 

Vi  v-'  / .  .  \  )  (1) 

w(ti,  tj) 

where  w(ti,tj)  is  a  weighting  function  that  monotonically  decreases  as  the  difference  between  ti 
and  tj  increases.  For  the  weighting  function,  we  use  the  following  definition  based  on  a  Gaussian 
kernel  with  standard  deviation  a, 

w{t,'ti)  =  7jSe:LV{J±wL)  •  <2) 

When  smoothing  horizontally,  we  use  a  =  5  s.  When  smoothing  vertically,  we  use  a  =  15  s.  A  larger 
a  is  required  for  vertical  smoothing  because  of  100  ft  Mode  C  quantization.  We  chose  these  values 
for  <7  after  trying  different  standard  deviations  on  a  sampling  on  horizontal  and  vertical  profiles  in 
our  data  set;  the  chosen  values  preserve  the  underlying  tracks  while  removing  noise. 

3.4  INTERPOLATION 

The  time  interval  between  radar  scans  in  our  data  is  much  longer  than  the  1  s  time  step  of 
our  dynamic  Bayesian  network.  Terminal  (short  range)  radars  scan  aircraft  approximately  every  5 
seconds,  and  en  route  (long  range)  radars  scan  aircraft  every  10  to  12  seconds.  Additionally,  it  is 
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common  for  sensors  to  skip  one  or  more  consecut  ive  scans  of  a  target  and  some  scans  produce  out  liers 
that  we  remove  (Section  3.2).  Hence,  we  need  to  use  interpolation  to  estimate  the  parameters  in 
our  dynamic  Bayesian  network.  We  chose  a  piecewise-cubic  Hermite  interpolation  scheme  that 
preserves  monotonicity  and  shape  [15]. 

Figure  6  shows  the  result  of  outlier  detection,  smoothing,  and  interpolation  on  a  track  from 
our  data  set.  Figure  7  shows  the  result  of  piecewise-cubic  Hermite  interpolation  on  a  smoothed 
track. 


3.5  FEATURE  EXTRACTION 

Feature  extraction  involves  converting  an  interpolated  track  into  sequences  of  airspace  classes, 
altitude  layers,  airspeeds,  vertical  rates,  turn  rates,  and  accelerations. 

•  Airspace  class:  We  estimate  latitude  and  longitude  of  radar  returns  using  an  algorithm 
developed  at  Lincoln  [16].  From  altitude  estimates  and  latitude  and  longitude  estimates,  we 
determine  the  class  of  airspace  by  searching  through  the  National  Airspace  System  Resources 
(NASR)  database  provided  by  the  FAA.  Since  the  altitude  estimates  are  based  on  Mode 
C  reports  of  pressure  altitude,  uncorrected  for  barometric  variation,  it  is  possible  that  the 
airspace  of  some  tracks  are  identified  incorrectly.  We  expect  that  some  inaccuracy  in  airspace 
class  identification  due  to  barometric  variation  has  a  negligible  impact  on  our  model. 

•  Altitude  layer:  Altitude  above  ground  level  (AGL)  determines  the  altitude  layer  in  our 
model.  We  estimate  altitude  AGL  by  subtracting  an  estimate  of  ground  elevation  from 
pressure  altitude.  Our  estimates  of  ground  elevation  come  from  Digital  Terrain  Elevation  Data 
(DTED)  provided  by  the  National  Geospatial-Intelligence  Agency  (NGA).  We  use  DTED 
Level  0,  which  has  post  spacing  of  30  arcseconds  (approximately  900  meters). 

•  Airspeed:  The  true  airspeed  at  time  t  is  given  by 

v(t)  =  \J (x(t  +  1)  -  i(i))2  4-  (y{t  +  1)  -  y(*))2  +  (h{t  +  1)  -  h{t))2  . 

•  Vertical  rate:  The  vertical  rate  is  estimated  from  the  smoothed  and  interpolated  altitudes 
estimated  from  Mode  C  reports.  The  vertical  rate  at  time  t  is  given  by  h(t)  =  h(t  +  1)  —  h(t). 

•  Turn  rate:  We  first  compute  the  heading  along  the  interpolated  track.  The  heading  at  time 
t  is  given  by  and  corresponds  to  the  direction  from  (x(t),  y(t))  to  (x(t  +  1  ),y{t  +  1)). 
To  compute  the  turn  rate  at  time  £,  we  find  the  acute  change  in  heading  between  ip(t)  and 
^{t+  1).  Turns  to  the  right  have  positive  turn  rates,  and  turns  to  the  left  have  negative  turn 
rates. 

•  Acceleration:  To  find  the  acceleration  at  a  particular  point,  we  average  the  change  in 
airspeed  per  unit  time  looking  forward  one  time  step  and  looking  back  one  time  step. 

Figure  8  shows  the  result  of  feature  extraction  on  the  same  track  shown  in  Figure  6. 
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TABLE  2.  Cut  points  used  for  feature  quantization. 


Cut  Points 

A 

B,  C,  D,  O 

L 

1200,3000, 5000 

V 

30,60,90,120,140,165,250 

V 

-1,-0.25,0.25,1 

h 

-1250,  -750,  -250, 250,  750, 1250 

ip 

-6,  -4.5, -1.5, 1.5, 4.5, 6 

3.6  FEATURE  SMOOTHING 

We  then  smooth  the  extracted  features  using  the  same  smoothing  scheme  we  used  for  tracks 
(Section  3.3).  For  turn  rate,  airspeed,  and  acceleration,  we  set  a  to  10  s,  20  s,  and  20  s  respectively. 
We  choose  these  numbers  large  enough  so  that  noise  is  removed  from  the  measurements  but  low 
enough  so  that  the  underlying  properties  of  the  maneuvers  are  not  lost.  We  do  not  smooth  vertical 
rates  in  this  step  because  the  altitudes  are  already  smoothed  (Section  3.3). 

3.7  QUANTIZATION 

In  order  to  be  modeled  by  a  discrete  Bayesian  network,  it  is  necessary  to  quantize  the  features. 
We  quantize  continuous  values  by  defining  a  sequence  of  cut  points  ci, . . .  ,cn.  Values  less  than  c\ 
are  in  the  first  bin,  values  greater  than  cn  are  in  the  (n  +  l)th  bin,  and  values  in  the  half-open 
interval  [cj_i,q)  are  in  the  ith  bin.  The  cut  points  we  used  for  quantization  are  listed  in  Table  2. 
The  cut  points  were  chosen  to  capture  the  variation  of  the  features  as  shown  in  the  histograms  in 
Figure  9. 

3.8  STATISTICS  EXTRACTION 

With  structures  for  the  initial  and  transition  distributions  and  the  quantized  features  from 
a  set  of  tracks,  we  are  able  to  collect  the  sufficient  statistics  to  estimate  the  parameters  for  our 
model.  For  the  two  Bayesian  networks,  the  sufficient  statistics  are  simply  the  counts  of  the 
various  features  (see  Appendix  C).2  Appendix  A  describes  the  sufficient  statistics  extracted  from 
our  data. 

Figure  10  illustrates  the  convergence  of  the  Bayesian  network  parameters.  The  horizontal  axis 
represents  the  number  of  samples  used  to  estimate  the  parameters  of  the  Bayesian  network,  and 
the  vertical  axis  represents  the  maximum  difference  between  elements  in  the  conditional  probability 
tables  with  the  addition  of  more  data.  As  the  figure  shows,  the  change  in  the  probabilities  in  the 
conditional  probabilities  converges  to  zero. 

2The  counts  are  called  sufficient  statistics  because  they  provide  a  summarization  of  the  data  that  is  sufficient  to 
compute  the  posterior  distribution  from  the  prior.  For  an  introduction  to  Bayesian  statistics,  see  [17]. 
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This  section  has  described  how  to  construct  a  model  of  nominal  VFR  flight  based  on  radar 
data.  The  next  section  describes  how  to  use  this  model  to  produce  new  trajectories  through 
sampling. 
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FIGURE  6.  Preprocessing.  Blue  lines  show  an  example  raw  track.  Red  lines  show  the  track  after  outlier  removal, 
smoothing,  and  interpolation. 
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FIGURE  7.  Piecewise-cubic  Hermite  interpolation  on  a  smoothed  track.  Red  crosses  indicate  smoothed  data 
points,  and  the  blue  curve  shows  the  interpolation. 
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FIGURE  8.  A  plot  of  extracted  features  over  time.  Blue  lines  show  features  before  smoothing,  and  red  lines 
show  features  after  smoothing.  The  green  blocks  show  the  outlines  of  the  bins  to  which  the  features  belong. 
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FIGURE  9.  Feature  histograms  of  observed  radar  data  based  on  97  million  samples. 
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Number  of  samples 

(a).  Initial  network. 


Number  of  samples 

(b).  Transition  network 


FIGURE  10.  Convergence  curve  for  Bayesian  network  conditional  probabilities. 
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4.  SAMPLING 


Once  we  process  the  data  as  described  in  the  previous  section,  we  can  use  the  model  structure 
and  sufficient  statistics  to  produce  new  tracks  that  are  representative  of  the  ones  observed  by  radar. 
The  first  step  involves  sampling  from  the  discrete  Bayesian  networks  representing  the  initial  and 
transition  distributions.  The  second  step  involves  converting  the  discrete  sample  into  a  continuous 
sample  by  sampling  within  bins.  Figure  11  illustrates  this  process.  An  example  of  working  through 
this  process  using  the  data  provided  for  this  encounter  model  is  given  in  Appendix  B. 


FIGURE  11.  Sampling  process  flow. 


4.1  DISCRETE  SAMPLING 

We  begin  by  sampling  from  the  Bayesian  network  representing  the  initial  distribution.  As 
outlined  in  Appendix  Cand  shown  in  more  detail  in  Appendix  B,  we  randomly  assign  value  k  to 
variable  X\  with  probability 

CXijk  h  Nijh 

Sfc'=i(ayfc'  Nijk') 

where 

•  j  is  the  instantiation  of  the  parents  of  Xj  in  the  Bayesian  network, 

•  Nijk  is  the  number  of  times  Xi  was  equal  to  k  when  its  parents  were  instantiated  to  j  in  the 
data, 
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•  otijk  is  a  Dirichlet  prior  parameter,  and 

•  Vi  is  the  number  of  ways  to  instantiate  the  parents  of  X{. 

We  use  OLijk  =  1,  which  corresponds  to  the  prior  assumption  that  all  combinations  of  relative 
frequencies  for  k  are  equally  probable.  Sampling  from  the  posterior  distribution  with  =  1  is 
equivalent  to  adding  1  to  all  the  counts  in  the  tables  in  Appendix  A  and  sampling  according  to  the 
resulting  relative  frequencies.  This  ensures  that  there  are  no  transitions  with  zero  probability  in 
our  Markov  model. 

Our  sample  from  the  Bayesian  network  tells  us  into  which  bins  airspeed,  vertical  rate,  turn 
rate,  and  acceleration  fall.  We  then  sample  within  the  bins  as  discussed  in  Section  4.2. 

Once  we  have  the  initial  state  of  the  trajectory,  we  can  sample  from  our  dynamic  Bayesian 
network  representing  how  the  state  changes.  We  fix  h(t),  and  v(t ),  and  then  use  the  standard 
Bayesian  network  sampling  scheme  to  determine  h(t  +  1),  ip(t  +  1),  and  v(t  4-  1).  The  process  may 
be  repeated  for  as  long  as  we  wish  to  run  the  trajectory. 


4.2  CONTINUOUS  SAMPLING 

To  produce  a  continuous  sample  from  the  discrete  sample  from  the  initial  distribution,  we 
simply  sample  uniformly  within  the  bins.  For  example,  if  we  determine  that  the  initial  airspeed 
is  within  the  bin  [60,90),  i.e.  the  second  bin  according  to  the  quantization  scheme  in  Table  2,  we 
simply  sample  from  the  uniform  distribution  over  the  half-open  interval  [60,90).  Because  the  first 
and  last  bins  associated  with  each  interval  are  unbounded  and  uniform  sampling  on  an  unbounded 
interval  is  undefined,  it  is  necessary  to  impose  some  bounds.  Table  3  shows  the  boundaries  in  our 
quantization  based  on  the  limits  observed  in  the  radar  data. 

TABLE  3.  Sampling  boundaries. 


Boundaries 

h 

500, 1200, 3000,  5000, 12500 

V 

0, 30, 60 

,90,120, 140,165,250, 300 

V 

-2,-1, 

-0.25,0.25,1,2 

h 

-2000,  - 

-1250,  -750,  -250, 250,  750, 1250, 2000 

-8,  -6, 

-4.5, -1.5, 1.5, 4.5, 6, 8 

To  produce  a  continuous  sample  of  the  trajectory  produced  by  the  transition  network,  we 
iterate  through  the  arrays  of  climb  rates,  turn  rates,  and  accelerations.  Instead  of  sampling  within 
bins  at  every  time  step,  we  only  produce  a  new  continuous  sample  at  some  fixed  rate  estimated 
from  the  data.  The  rates  we  used  are  0.0211495,  0.0412835,  and  0.0885336  changes  per  second  for 
acceleration,  vertical  rate,  and  turn  rate  respectively.  We  estimated  these  rates  from  the  data  by 
introducing  3  smaller  bins  within  each  bin  and  computing  the  relative  frequency  that  tracks  stay 
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within  a  single  smaller  bin  (versus  moving  to  another  small  bin  within  the  same  coarse  bin).  A 
similar  strategy  was  used  by  Eurocontrol  for  their  cooperative  encounter  model. 

When  producing  continuous  samples  from  bins  that  include  zero  in  their  range,  we  simply 
produce  zero  instead  of  sampling  uniformly  in  order  to  prevent  very  small  vertical  rates,  turn  rates, 
and  acceleration. 

Figure  12  plots  an  example  of  vertical  rates,  turn  rates,  accelerations,  and  airspeeds  generated 
by  sampling  from  the  Bayesian  networks.  Figure  13  shows  the  resulting  track  produced  by  the 
sampled  features  shown  in  the  previous  figure.  Translation  of  features  into  tracks  involves  running 
a  discrete-time  simulation,  described  in  the  next  section. 


0.5 


a> 

o 

u 

< 


-0.5 


200  400 

Time  (s) 


600 


FIGURE  12.  A  plot  of  sampled  features  over  time.  Blue  lines  show  continuous  samples,  and  green  blocks  show 
the  outlines  of  the  discrete  samples. 
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FIGURE  13.  A  track  generated  by  sampling  from  the  initial  and  transition  distributions.  Positions  and  altitudes 
are  relative  to  the  initial  position  and  altitude. 
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5.  SIMULATION 


Earlier  sections  explained  how  to  build  a  dynamic  probabilistic  model  of  aircraft  and  how  to 
sample  from  the  model  to  produce  nominal  trajectories  that  are  representative  of  what  we  observed 
in  our  radar  data.  The  same  modeling  techniques  discussed  earlier  may  be  used  to  build  a  model 
of  the  manned  or  unmanned  aircraft  with  the  collision  avoidance  system  to  be  evaluated.  This 
section  explains  how  to  combine  a  model  of  the  aircraft  with  the  collision  avoidance  system  to  be 
evaluated,  which  we  call  AC1,  with  a  model  of  an  intruder,  which  we  call  AC2,  into  a  simulated 
encounter. 

The  trajectory  for  AC1  may  be  specified  by  the  analyst  (e.g.,  to  focus  on  a  particular  phase 
of  flight),  based  on  actual  flight  paths  from  mission  planning  or  radar  data,  or  randomly  generated 
using  a  statistical  model  representative  of  that  aircraft’s  typical  flight  profiles.  In  a  study  of  VFR- 
011-VFR  encounters,  for  example,  trajectories  for  both  AC1  and  AC2  could  be  generated  using  the 
uncorrelated  encounter  model  described  here. 

We  avoid  simulating  encounters  that  are  extremely  unlikely  to  result  in  NMACs  by  focusing 
our  computational  effort  on  encounters  that  occur  in  an  encounter  cylinder  centered  on  AC1.  AC2 
is  initialized  on  the  surface  of  the  encounter  cylinder  and  the  dynamic  models  are  used  to  update 
the  states  of  AC1  and  AC2  over  time.  If  the  collision  avoidance  system  on  AC1  or  AC2  suggests  an 
avoidance  maneuver,  it  overrides  the  nominal  rates  commanded  by  the  encounter  model.  If  AC2 
enters  an  NMAC  cylinder  or  exits  the  encounter  cylinder,  the  encounter  run  is  terminated.  The 
NMAC  cylinder  has  radius  rnmac  and  height  2/*nmac.  We  use  rnmac  =  500ft  and  hn mac  =  100  ft. 

Running  millions  of  randomly  generated  encounters  allows  us  to  estimate  the  probability 
that  an  aircraft  that  enters  the  encounter  cylinder  penetrates  the  NMAC  cylinder  before  exiting 
the  encounter  cylinder.  We  use  P(nmac  |  enc)  to  represent  this  probability.  Section  6  explains  how 
to  use  this  probability  to  estimate  the  mean  time  between  NMACs. 


5.1  ENCOUNTER  CYLINDER  DIMENSIONS 

The  encounter  cylinder  lias  radius  renc  and  height  2 henc.  The  appropriate  dimensions  of  the 
encounter  cylinder  depend  on  the  aircraft  dynamics  and  collision  avoidance  system.  If  the  encounter 
cylinder  is  too  small,  the  collision  avoidance  system  will  not  have  enough  opportunity  to  be  fully 
exercised  before  a  collision.  If  the  encounter  cylinder  is  too  large,  then  computation  is  wasted. 

An  upper  bound  for  renc  is  the  amount  of  time  required  by  the  collision  avoidance  system  to 
detect  and  track  a  target  multiplied  by  the  sum  of  the  maximal  airspeeds  of  AC1  and  AC2.  An 
upper  bound  for  henc  is  the  amount  of  time  required  by  the  collision  avoidance  system  to  detect 
and  track  a  target  multiplied  by  the  sum  of  the  maximal  vertical  rate  magnitudes  of  ACl  and  AC2. 
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hTe]  <  0 


hTe'  >  0 


R  ■  (v2  -  v\)  <  0 


FIGURE  14.  Depending  on  where  on  the  encounter  cylinder  the  intruder  is  initialized,  different  criteria  must  be 
met  in  order  to  accept  a  given  trajectory  for  further  simulation. 


5.2  INITIAL  CONDITIONS 

We  use  rejection  sampling  to  generate  the  initial  conditions  of  an  encounter.  Rejection  sam¬ 
pling  involves  proposing  a  series  of  candidate  samples  from  a  random  distribution  until  choosing 
one  that  meets  a  set  of  criteria  (summarized  in  Figure  14).  The  process  we  use  for  generating  initial 
conditions  for  encounters  is  as  follows: 

1.  Generate  airspeeds,  vertical  rates,  turn  rates,  and  accelerations  for  AC1  and  AC2  according 
to  their  models. 

2.  Randomly  initialize  the  position  of  AC2  on  the  surface  of  the  encounter  cylinder  centered  on 
ACl.  AC2  is  not  restricted  to  being  on  the  side  of  the  cylinder;  it  may  be  initialized  on  the 
top  or  bottom.  The  probability  of  being  initialized  on  the  top,  bottom,  or  side  is  proportional 
to  its  surface  area.  The  bearing  of  AC2  relative  to  ACl  is  denoted 

3.  The  heading  of  ACl  is  set  to  0.  The  heading  of  AC2,  denoted  0,  is  randomly  selected  from 
a  uniform  distribution  over  [0,27r). 

4.  If  AC2  was  initialized  on  the  top  of  the  encounter  cylinder,  accept  the  sample  if  the  vertical 
rate  of  AC 2  relative  to  AC  1,  denoted  hTe\  is  negative.  This  ensures  that  AC2  is  penetrating 
the  encounter  cylinder  in  the  first  time  step. 

5.  If  AC2  was  initialized  on  the  bottom  of  the  encounter  cylinder,  accept  the  sample  if  the 
vertical  rate  of  AC 2  relative  to  AC  1,  denoted  hre\  is  positive.  This  ensures  that  AC2  is 
penetrating  the  encounter  cylinder  in  the  first  time  step. 

6.  If  AC2  was  initialized  on  the  side  of  the  encounter  cylinder,  accept  the  sample  if  R  •  (V2  —  v  i ) 
is  negative,  where  R  =  (sin  x,  cos  x)  is  the  unit  vector  from  ACl  to  AC2.  The  vectors  V\ 
and  V2  are  the  velocities  of  ACl  and  AC2  respectively.  When  R  •  (v2  —  V\)  is  negative,  the 
relative  velocity  of  AC2  is  into  the  encounter  cylinder,  and  therefore  should  be  accepted.  See 
Figure  15. 
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V2  ~  Vi 


FIGURE  15.  Variables  relevant  to  rejection  sampling  in  the  horizontal  plane. 

The  process  is  repeated  until  a  candidate  initialization  is  accepted. 

A  byproduct  of  rejection  sampling  is  that  the  bearing  distribution  is  nonuniform.  As  one  would 
expect,  more  encounters  occur  head  on  than  from  the  side.  Figure  16  shows  bearing  distributions 
of  VFR/VFR  encounters. 

5.3  TRAJECTORY  CONSTRUCTION 

Once  the  initial  conditions  are  selected,  the  dynamic  models  of  AC1  and  AC2  are  used  to 
update  their  trajectories  over  time. 

Given  the  initial  state  of  an  aircraft  (including  at  a  minimum:  position,  altitude,  heading, 
climb  rate,  turn  rate,  and  velocity)  and  its  control  variables  (h,  ip,  and  v)  the  aircraft  state  is 
updated  in  1  s  time  steps  using  a  dynamic  model.  Due  to  the  wide  variety  of  possible  dynamic 
model  implementations,  details  for  computing  this  state  update  are  not  provided  here.  A  basic 
approach  would  be  to  apply  simple  point-mass  kinematics  to  update  the  aircraft  state  without 
considering  what  the  aircraft  may  actually  be  doing  in  terms  of  bank  angle,  pitch  rate,  etc.  More 
complex  implementations  could  include  6  degree-of- freedom  dynamic  models  in  which  the  control 
variables  are  treated  as  target  states  provided  to  an  autoflight  control  system  which  then  applies 
the  necessary  control  deflections  to  attain  those  targets.  Lincoln  Laboratory’s  Collision  Avoidance 
System  Safety  Assessment  Tool  (CASSATT)  typically  uses  a  4  degree-of-freedom  model  to  update 
aircraft  state  by  applying  the  necessary  airspeed  acceleration,  roll  rate,  and  pitch  rate  to  achieve 
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FIGURE  16.  Initial  bearing  distributions  produced  by  rejection  sampling  for  two  randomly-generated  VFR 
aircraft. 

the  target  values  for  h,  and  v  assuming  curvilinear  motion  with  a  zero-sideslip  constraint.  The 
zero-sideslip  constraint  can  be  relaxed,  for  instance,  if  it  is  necessary  to  model  transient  dynamics 
with  higher  fidelity. 

A  simulation  run  terminates  when  the  intruder  either  exits  the  encounter  cylinder  or  pene¬ 
trates  the  NMAC  cylinder.  After  running  many  simulations,  we  estimate  P(nmac  |  enc)  by  dividing 
the  number  of  runs  that  resulted  in  an  NMAC  by  the  total  number  of  runs. 


5.4  MULTIPLE  ENCOUNTERS 


The  approach  discussed  in  Section  5.2  only  handles  pairwise  encounters.  The  probability 
of  an  intruder  penetrating  the  encounter  cylinder  while  another  intruder  is  within  the  encounter 
cylinder  is  likely  to  be  very  small.  One  may  compute  this  probability  using 


(4) 


where  p(t)  is  the  distribution  over  the  amount  of  time  intruders  spend  in  the  encounter  cylinder  and 
Aenc  is  the  rate  at  which  new  intruders  penetrate  the  encounter  cylinder.  In  the  above  equation, 
l_eAenc*  COmes  from  the  cumulative  density  function  for  an  exponential  distribution.  The  possibility 
of  simultaneous  multiple  intruders  needs  to  be  examined  and  will  be  an  area  of  future  work. 


28 


6.  SAFETY  EVALUATION 


The  previous  section  described  how  to  compute  P(nmac  |  enc),  the  probability  that  an  aircraft 
that  enters  the  encounter  cylinder  penetrates  the  NMAC  cylinder  before  exiting  the  encounter 
cylinder.  This  section  explains  how  to  compute  the  rate  at  which  aircraft  enter  the  encounter 
cylinder,  Aenc-  Once  we  know  P(nmac  |  enc)  and  Aenci  we  may  estimate  the  NMAC  rate 


Anmac  —  P(nmac  |  enc)Aenc  . 

The  mean  time  between  NMACs  is  simply  Am]iac. 

We  begin  with  the  assumption  that  the  density  of  traffic  outside  the  encounter  cylinder  is 
fixed  and  known.  This  density,  p,  has  units  (NM2  x  ft)-1.  The  rate  at  which  aircraft  enter  the 
encounter  cylinder  depends  on  p  and  the  average  volume  of  new  airspace  the  encounter  cylinder 
sweeps  through  per  unit  time. 

If  the  average  relative  horizontal  speed  is  sj^Jriz  and  the  average  relative  vertical  speed  is  s^rt, 
then  the  average  volume  of  new  airspace  the  encounter  cylinder  sweeps  through  per  unit  time  is 

‘Irenc/iencShcH-iz;  ^*^enc^vert  *  ('*) 


The  average  relative  vertical  speed  is  given  by 


srel 

^vert 


/oo  re* 

-oc  J  —c 


p(hi)p(h2) 


hi  —  h2 


d/?  i  d/?2  5 


(6) 


where  h\  is  the  vertical  rate  of  the  unmanned  aircraft  and  /i2  is  the  vertical  rate  of  the  intruder. 
The  average  relative  horizontal  speed  is  given  by 


-rel 

'Uioriz 


i  r  7r  roc  roc  , - 

=  —  J  J  J  p(vf  )p(vg )  y  -I-  cos  a)2  +  sin  a)2  dt^dt’fda  , 


(7) 


where  v\  is  the  groundspeed  of  the  unmanned  aircraft  and  vg2  is  the  ground  speed  of  the  intruder. 
Both  Equations  6  and  7  may  be  estimated  through  sampling.  We  estimate  5{^Jriz  to  be  147  kt  for 
VFR/VFR  encounters.  We  estimate  sT^rt  to  be  260ft/min  for  VFR/VFR  encounters. 

If  we  make  the  following  two  assumptions: 


1.  the  density  of  air  traffic  outside  the  encounter  cylinder  is  uniform,  and 

2.  the  trajectories  of  aircraft  outside  of  the  encounter  cylinder  are  independent  of  the  trajectories 
of  aircraft  within  the  encounter  cylinder, 

then  the  encounter  rate  is 

Aenc  =  P  ^drenc^enc^horiz  +  ^enc^vert)  •  (8) 

The  NMAC  rate  is  simply  the  product  of  Aenc  and  P(nmac  |  enc). 
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7.  SUMMARY  AND  FURTHER  WORK 


This  report  has  presented  a  new  approach  to  encounter  modeling  that  allows  for  the  generation 
of  more  realistic  encounters  than  previous  models.  The  approach  involves  modeling  the  dynamics 
of  aircraft  state  based  on  a  Markov  model  where  the  probability  of  the  next  state  depends  only 
upon  the  current  state.  One  way  to  represent  a  Markov  model  is  with  an  exhaustive  state-transition 
matrix  that  specifies  the  probability  of  transitioning  between  all  pairs  of  states.  Unfortunately,  the 
number  of  independent  parameters  required  to  define  the  matrix  grows  super  exponentially  with 
the  number  of  variables  defining  the  model.  The  more  independent  parameters  there  are  in  the 
model,  the  more  data  one  needs  to  properly  estimate  their  values.  However,  by  using  dynamic 
Bayesian  networks,  we  can  leverage  conditional  independence  between  some  variables  to  greatly 
reduce  the  number  of  parameters.  We  can  learn  the  structure  of  the  dynamic  Bayesian  network  by 
maximizing  the  posterior  probability  of  the  network  structure  given  the  data. 

The  work  presented  in  this  report  has  focused  on  a  model  where  the  trajectories  of  the 
aircraft  involved  in  the  encounter  are  independent  of  each  other  prior  to  intervention  by  a  collision 
avoidance  system,  human  or  automated.  This  model  assumes  that  aircraft  blunder  into  close 
proximity  without  prior  intervention.  We  are  currently  developing  a  separate  model  for  aircraft 
under  air  traffic  control,  where  intervention  may  impact  the  way  in  which  aircraft  encounter  each 
other.  Such  a  correlated  encounter  model  is  similar  to  the  uncorrelated  encounter  model  described 
here  except  that  the  variables  defining  the  state  of  the  aircraft  in  the  encounter  are  tied  to  each  other 
in  the  dynamic  Bayesian  network.  We  are  also  developing  a  noncooperative  model  of  unconventional 
aircraft  such  as  balloons,  ultralights,  and  gliders. 

The  uncorrelated  encounter  model  presented  in  this  report  along  with  a  correlated  encounter 
model  under  development  will  be  used  to  generate  encounter  situations  for  use  in  Monte  Carlo 
safety  analyses  of  collision  avoidance  systems  for  manned  and  unmanned  aircraft.  The  results  of 
these  robustness  studies  will  inform  the  development  and  certification  of  new  systems. 

The  model  presented  in  this  report  will  be  reevaluated  over  Spring  and  Summer  2008,  and  a 
final  version  will  be  published  in  September  2008. 
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APPENDIX  A 
MODEL  PARAMETERS 


This  appendix  describes  the  sufficient  statistics,  Nijk  (see  Appendix  C),  used  to  estimate  the 
conditional  probabilities  associated  with  the  initial  and  transition  distributions.  These  sufficient 
statistics  are  based  on  beacon  reports  from  aircraft  squawking  VFR  (Mode  A  code  1200)  between 
December  1,  2007  and  December  7,  2007.  This  amounts  to  30,109,411  reports  representing  78,000 
flight  hours  from  across  the  United  States.  Other  parameters  that  are  relevant  to  generating  new 
encounters  from  the  model  are  also  described  in  this  section. 

A  text  file,  available  electronically  from  Lincoln  Laboratory,  describes  the  following  model 
parameters: 

•  Variable  labels:  A  quoted,  comma-delimited  list  specifies  the  variable  labels,  e.g.  \dot  \psi, 
as  would  be  used  by  DT£X.  There  are  different  variable  labels  for  the  initial  network  and  the 
transition  network. 

•  Graphical  structure:  A  binary  matrix  is  used  to  represent  graphical  structure.  A  1  in  the 
zth  row  and  jth  column  means  that  there  is  a  directed  edge  from  the  ith  variable  to  the  jth 
variable  in  the  Bayesian  network.  The  text  file  specifies  two  graphical  structures;  one  for  the 
initial  network  and  the  other  for  the  transition  network.  The  element  in  the  ith  row  and  j th 
column  is  represented  by  G(i,  j). 

•  Variable  instantiations:  For  each  network,  a  list  of  integers  specifies  the  number  of  instan¬ 
tiations  that  exist  for  each  variable. 

•  Sufficient  statistics:  For  each  network,  a  list  of  integers  specifies  the  sufficient  statistics. 
We  explain  how  to  interpret  the  array  of  integers  below. 

•  Boundaries:  The  boundaries  of  the  variable  bins  are  specified  by  a  row  of  numbers.  The 
variables  A  and  L  are  not  quantized  because  they  are  already  discrete,  and  so  boundaries  do 
not  exist.  A  *  is  used  for  these  variables. 

•  Resampling  rates:  A  list  of  numbers  specify  the  resampling  rates  (Section  4.2). 

The  list  of  numbers  describing  the  sufficient  statistics,  A^*.,  requires  explanation.  The  array 
is  ordered  first  by  increasing  fc,  then  increasing  j,  and  then  increasing  i.  For  some  specified  variable 
parental  instantiation  j,  and  variable  instantiation  k,  the  value  Ay*  is  given  by  the  following 
element  in  the  list 

t-i 

k  +  n(j  -  1)  +  gyry  ,  (A-l) 

i'= 1 

where  q  and  r  are  as  specified  in  Appendix  C.  It  is  important  to  clarify  the  ordering  of  the  parental 
instantiations.  If  the  variables  are  instantiated  to  bins  the  parental  instantiation  of 
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variable  i  is  given  by 


n 


i'  —  l 


j  —  i  +  o(^*'  ~  i)  U 


(A-2) 


The  following  is  a  fragment  of  the  parameter  file.  The  lines  that  describe  the  sufficient 
statistics  are  truncated  due  to  length. 

#  labels.initial 

"A",  "L",  "v",  "\dot  v" ,  "\dot  h" ,  "\dot  \psi" 

#  G_initial 
0  11111 
0  0  1111 
0  0  0  1  1  1 
0  0  0  0  1  1 
0  0  0  0  0  1 
0  0  0  0  0  0 

#  r.initial 
4  4  8  5  7  7 

#  N.initial 

487885  517728  10285433  86339819  249426  130407  78279  29773  302812  182518  . . . 

#  labels.transition 

"A",  "L",  "v",  "\dot  v(t)u,  "\dot  h(t) " ,  "\dot  \psi(t)u, 

"\dot  v(t+l) " ,  "\dot  h(t+l) " ,  "\dot  \psi(t+l)" 

#  G_transition 
000000011 
0  0  0  0  0  0  1  1  1 
000000000 
000000100 
000000010 
000000001 
000000000 
0  0  0  0  0  0  1  0  1 
000000100 

#  r_transition 
448577577 

#  N.transition 

135  1  0  0  0  168  1  0  0  0  39  0  0  0  0  59  0  0  0  0  3  921  46  0  0  3  902  20  0  0  ... 


#  boundaries 


0  30  60  90  120  140  165  250  300 
-2-1-0012 

-2000  -1250  -750  -250  250  750  1250  2000 
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-8-6-5-22568 
#  resample.rates 

000  0.0211495  0.0412835  0.0885336 


This  document  and  the  parameters  file  is  available  for  download  from  ftp://duke.atc.ll. 
mit.edu  as  a  ZIP  file  located  at  /pub/outgoing/encounter_model/uncor-20080303 . zip.  This 
FTP  site  allows  for  anonymous  access.  Sign  in  with  the  username  anonymous  and  your  email 
address  as  your  password. 
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APPENDIX  B 

TRAJECTORY  GENERATION 

This  appendix  explains  how  to  generate  a  trajectory  from  our  model. 


B.l  INITIAL  NETWORK  SAMPLING 


The  first  step  in  generating  a  random  trajectory  using  our  model  is  to  sample  from  the 
Bayesian  network  representing  the  initial  state  distribution.  To  sample  from  a  Bayesian  network, 
as  explained  in  Appendix  C,  is  to  produce  a  topological  sort  of  the  nodes  in  the  network.  A 
topological  sort  orders  the  nodes  of  the  network  so  that  parents  come  before  their  descendants. 
The  following  is  the  graphical  structure  of  the  initial  network  as  specified  in  the  parameters  file 
(see  Appendix  A)  and  shown  in  Figure  B-l: 


0  11111 
0  0  1111 
0  0  0  1  1  1 
0  0  0  0  1  1 
0  0  0  0  0  1 
0  0  0  0  0  0 


As  can  be  seen,  this  ordering  of  the  nodes  is  already  topologically  sorted:  the  first  node 
(airspace  class)  is  connected  to  all  other  nodes.  The  second  node  (altitude  layer)  is  connected  to 
all  following  nodes,  and  so  on.  The  final  node  (0)  is  not  connected  to  any  other  nodes.  We  will  see 
later  that  the  transition  network,  in  contrast,  requires  a  topological  sort. 

With  the  nodes  being  topologically  sorted,  we  begin  by  sampling  the  first  variable.  As  specified 
in  the  parameters  file,  the  first  variable  is  A,  airspace  class.  Equation  C-5,  reproduced  below,  shows 
how  to  produce  a  random  sample: 


P(Xi  =  k\nij,D,G)  = 


&ijk  T  Nijk 
Ylk'=i(aijk'  T  Njjk') 


In  other  words,  the  probability  of  selecting  a  particular  airspace  class  is  proportional  to  the  prior 
(otijk)  plus  the  frequency  that  airspace  class  appeared  in  the  data  (iV^*).  Following  Cooper  and 
Herskovits  [18],  we  use  an  objective  prior  and  set  to  1.  To  determine  the  values  for  NtJ^,  we 
need  to  look  at  the  sufficient  statistics  recorded  in  the  parameters  file.  The  sufficient  statistics 
for  the  initial  network  is  recorded  as  a  long  series  of  numbers.  We  use  the  scheme  described  in 
Appendix  A  to  determine  the  actual  values.  These  values  turn  out  to  be  the  first  four  numbers  in 
the  sufficient  statistics  sequence  and  are  shown  in  the  following  table. 
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FIGURE  B-l.  A  graphical  representation  of  the  initial  network. 
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TABLE  B-l.  N(A) 


A 

B 

C 

D 

O 

487885 

517728 

10285433 

86339819 

Now,  we  may  compute  the  probability  of  selecting  the  different  values  of  A  using  Equation  C-5: 


•  P(A  =  B)  =  (487885  +  l)/(487885  +  1  +  517728  +  1  4-  10285433  +  1  +  86339819  +  1)  = 
487886/97630869  =  0.00499725143284344 

•  P{A  =  C)  =  (517728  4-  l)/(487885  +  1  4-  517728  4-  1  4-  10285433  +  1  4-  86339819  +  1)  = 
517729/97630869  =  0.00530292319737521 

•  P(A  =  D)  =  (10285433  4-  l)/(487885  +  1  4-  517728  4-  1  4-  10285433  4-  1  4-  86339819  4-  1)  = 
10285434/97630869  =  0.105350224835139 

•  P{A  =  O)  =  (86339819  4-  l)/(487885  4-  1  4-  517728  4-  1  4-  10285433  4-  1  4-  86339819  +  1)  = 
86339820/97630869  =  0.884349600534642 


We  use  a  random  number  generator  to  choose  an  airspace.  For  this  example,  let  us  choose  (),  the 
fourth  instantiation  of  A. 

The  next  step  is  to  instantiate  the  second  variable,  L,  which  is  altitude  layer.  Choosing  a 
random  instantiation  for  L  is  not  quite  as  easy  as  for  A  because  L  depends  upon  other  variables, 
namely  A.  We  need  to  compute  the  conditional  probability  distribution  P(L  \  A  =  O),  i.e.  the 
distribution  over  the  values  of  L  given  that  A  is  O.  We  consult  the  sufficient  statistics  table  N(L  \  A) 
which  is  extracted  from  the  parameters  file  (using  the  process  explained  in  Appendix  A)  and  is 
displayed  in  the  table  below. 


TABLE  B-2. 

N(L  |  A) 

A 

L 

[500,  1200) 

[1200,  3000) 

[3000.  5000) 

[5000,  oo] 

B 

249426 

130407 

78279 

29773 

C 

302812 

182518 

32339 

59 

D 

8383276 

1901858 

299 

0 

O 

27280449 

39197911 

13704707 

6156752 

Since  we  are  only  interested  in  the  counts  associated  with  A  =  O,  we  consider  only  the  last 
row  in  the  table.  To  translate  the  counts  into  probabilities,  we  add  1  to  each  element  in  the  last  row 
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and  then  divide  by  the  total  resulting  sum  of  that  row.  Again,  we  use  a  random  number  generator 
to  select  an  airspace  layer  according  to  these  probabilities.  This  process  of  randomly  assigning 
values  to  each  of  the  variables  conditional  on  the  values  of  their  parents  continues  until  all  of  the 
variables  have  been  assigned.  For  this  example,  assume  that  the  following  assignments  have  been 
made: 

•  A  =  0 

•  L  =  [3000,  5000) 

•  v  =  [120, 140) 

•  v  =  [-0.25,0.25) 

•  h  =  [250,  750) 

•  ip  =  [-1.5, 1.5) 


Specific  values  within  each  bin  are  then  determined  based  on  a  different  uniform  random 
number  for  each.  In  the  case  where  a  bin  spans  0  (as  it  does  in  this  example  for  u,  and  ^),  the 
value  0  itself  is  assigned. 


B.2  TRANSITION  NETWORK  SAMPLING 


We  have  described  how  to  sample  from  the  initial  Bayesian  network  to  generate  a  random  ini¬ 
tial  state.  The  next  step  is  to  use  the  dynamic  Bayesian  network  representing  the  transition  distribu¬ 
tion  to  generate  the  state  at  the  next  time  step  based  on  the  initial  state.  The  parameters  file  defines 
the  following  sequence  of  variables  in  the  dynamic  Bayesian  network:  A,  L,  t;,  v(t ),  h(t),  ^(t),  v(t  + 
l),/i(£  +  1  )1ip(t  +  1).  The  parameters  file  also  defines  the  following  graphical  structure  for  the 
network: 


000000011 

000000111 

000000000 

000000100 

000000010 

000000001 

000000000 

000000101 

000000100 


The  graphical  representation  of  this  matrix  is  shown  in  Figure  B-2.  As  can  be  seen,  the 
ordering  in  the  parameters  file  is  not  topologically  sorted  because  there  are  arrows  from  h(t  +  1) 
and  ^(t  +  1)  to  v(t  +  1)  but  v(t  +  1)  comes  before  both  h(t  +  1)  and  ip(t  -hi)  in  the  list.  We  need  to 
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FIGURE  B-2.  A  graphical  representation  of  the  transition  network. 


topologically  sort  the  list  of  variables  to  make  sure  that  we  assign  values  to  variables  in  the  right 
order.  One  possible  topological  sort  involves  simply  moving  v(t  4-  1)  to  the  end  of  the  list. 

We  are  only  interested  in  assigning  new  values  to  the  dynamic  variables,  namely  h(t  +  1), 
ip(t  +  1)  and  v(t+  1).  The  process  is  similar  to  the  process  used  to  sample  from  the  initial  network. 
First,  we  consult  the  table  of  sufficient  statistics  for  h(t  -I-  1).  This  table  of  conditional  counts  is 
shown  below.  For  sampling  purposes,  we  are  only  interested  in  the  row  that  represents  the  current 
variable  assignment,  i.e.  A  =  O,  L  =  [3000,  5000),  and  h  =  [250,  750).  The  relevant  line  highlighted 
in  the  table: 


TABLE  B-3.  N(h(t  +  1)  |  A,  L,  h(t)) 


h(t  +  1) 


A 

L 

h(t) 

[-2000,  -1250) 

[-1250,  -750) 

[-750,  -250)  [-250, 

250) 

[250,  750) 

[750.  1250) 

[1250,  2000] 

B 

[500,  1200) 

[-2000, 

-1250) 

437 

13 

0 

0 

0 

0 

0 

C 

[500,  1200) 

[-2000, 

-1250) 

1436 

23 

0 

0 

0 

0 

0 

D 

[500,  1200) 

[-2000. 

-1250) 

32319 

1059 

0 

0 

0 

0 

0 

O 

[500,  1200) 

[-2000. 

-1250) 

44854 

1503 

2 

0 

0 

0 

0 

B 

[1200,  3000) 

[-2000, 

-1250) 

576 

14 

0 

0 

0 

0 

0 

C 

[1200,  3000) 

[-2000, 

-1250) 

1202 

45 

0 

0 

0 

0 

0 

D 

[1200,  3000) 

(-2000, 

-1250) 

14412 

450 

0 

0 

0 

0 

0 

O 

[1200,  3000) 

[-2000, 

-1250) 

73685 

2484 

2 

0 

0 

0 

0 

B 

[3000,  5000) 

[-2000. 

-1250) 

1137 

34 

0 

0 

0 

0 

0 

C 

[3000,  5000) 

[-2000. 

-1250) 

551 

16 

0 

0 

0 

0 

0 

D 

[3000,  5000) 

[-2000. 

-1250) 

32 

1 

0 

0 

0 

0 

0 

O 

[3000,  5000) 

[-2000, 

-1250) 

49962 

1567 

2 

0 

0 

0 

0 

B 

(5000,  oo) 

[-2000. 

-1250) 

1452 

16 

0 

0 

0 

0 

0 

C 

[5000,  oo] 

[-2000, 

-1250) 

0 

0 

0 

0 

0 

0 

0 

D 

[5000,  oo] 

[-2000. 

-1250) 

0 

0 

0 

0 

0 

0 

0 

O 

[5000,  oo] 

[-2000, 

-1250) 

75599 

1325 

0 

0 

0 

0 

0 

B 

[500,  1200) 

[-1250, 

-750) 

23 

2229 

66 

0 

0 

0 

0 

C 

[500,  1200) 

[-1250, 

-750) 

48 

3693 

81 

0 

0 

0 

0 

D 

(500,  1200) 

[-1250, 

-750) 

1318 

176114 

5157 

0 

0 

0 

0 

O 

[500,  1200) 

[-1250, 

-750) 

1807 

259491 

8025 

0 

0 

0 

0 

B 

[1200,  3000) 

[-1250, 

-750) 

18 

1490 

55 

0 

0 

0 

0 

C 

[1200,  3000) 

[-1250, 

-750) 

52 

5107 

187 

0 

0 

0 

0 

Continued  on  next  page.  . 
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TABLE  B-3.  Continued 


A 

L 

hit) 

[-2000,  -1250) 

[-1250,  -750) 

[-750,  -250) 

[-250,  250) 

[250,  750) 

[750,  1250) 

[1250.  2000] 

D 

[1200,  3000) 

[-1250,  -750) 

494 

58181 

1782 

0 

0 

0 

0 

O 

[1200,  3000) 

[-1250,  -750) 

2726 

314132 

10417 

1 

0 

0 

0 

B 

[3000,  5000) 

[-1250,  -750) 

42 

1445 

82 

0 

0 

0 

0 

C 

[3000,  5000) 

[-1250,  -750) 

15 

865 

39 

0 

0 

0 

0 

D 

[3000,  5000) 

[-1250,  -750) 

1 

0 

1 

0 

0 

0 

0 

O 

[3000,  5000) 

[-1250,  -750) 

1620 

167085 

4740 

0 

0 

0 

0 

B 

[5000,  oo] 

[-1250,  -750) 

22 

428 

16 

0 

0 

0 

0 

C 

[5000,  oo] 

[-1250,  -750) 

0 

0 

0 

0 

0 

0 

0 

D 

[5000,  oo] 

[-1250,  -750) 

0 

0 

0 

0 

0 

0 

0 

O 

[5000,  oo] 

[-1250,  -750) 

1345 

111865 

2811 

1 

0 

0 

0 

B 

[500,  1200) 

[-750,  -250) 

0 

100 

15290 

186 

0 

0 

0 

C 

[500,  1200) 

[-750,  -250) 

0 

160 

37490 

500 

0 

0 

0 

D 

[500,  1200) 

[-750,  -250) 

0 

7037 

1712037 

23271 

0 

0 

0 

O 

[500,  1200) 

[-750,  -250) 

2 

10383 

2919815 

41529 

0 

0 

0 

B 

[1200,  3000) 

[-750,  -250) 

0 

61 

11823 

217 

0 

0 

0 

C 

[1200,  3000) 

[-750,  -250) 

0 

209 

45446 

639 

0 

0 

0 

D 

[1200,  3000) 

[-750,  -250) 

2 

1997 

388905 

5752 

0 

0 

0 

O 

[1200,  3000) 

[-750,  -250) 

1 

11662 

2796178 

41758 

0 

0 

0 

B 

[3000,  5000) 

[-750,  -250) 

0 

91 

9007 

165 

0 

0 

0 

C 

[3000,  5000) 

[-750,  -250) 

0 

41 

4714 

103 

0 

0 

0 

D 

[3000,  5000) 

[-750,  -250) 

0 

1 

30 

1 

0 

0 

0 

O 

[3000,  5000) 

[-750,  -250) 

1 

5082 

1001108 

12656 

0 

0 

0 

B 

[5000,  oo] 

[-750,  -250) 

0 

23 

1193 

28 

0 

0 

0 

C 

[5000,  oo] 

[-750,  -250) 

0 

0 

0 

0 

0 

0 

0 

D 

[5000,  oo] 

[-750,  -250) 

0 

0 

0 

0 

0 

0 

0 

O 

[5000,  oo] 

[-750,  -250) 

2 

2937 

490952 

4873 

0 

0 

0 

B 

[500,  1200) 

[-250,  250) 

0 

0 

380 

128943 

265 

0 

0 

C 

[500,  1200) 

[-250,  250) 

0 

0 

1021 

300731 

632 

0 

0 

D 

[500,  1200) 

[-250,  250) 

0 

0 

39831 

11486322 

20842 

0 

0 

O 

[500,  1200) 

[-250,  250) 

0 

0 

66410 

22520261 

36293 

1 

0 

B 

[1200,  3000) 

[-250,  250) 

0 

0 

247 

79669 

179 

0 

0 

C 

[1200,  3000) 

[-250,  250) 

0 

0 

840 

267771 

482 

0 

0 

D 

[1200,  3000) 

[-250,  250) 

0 

0 

7390 

2275618 

4651 

0 

0 

O 

[1200,  3000) 

[-250,  250) 

0 

0 

52430 

18977463 

29588 

1 

0 

B 

[3000.  5000) 

[-250,  250) 

0 

0 

187 

61715 

161 

0 

0 

C 

[3000,  5000) 

[-250,  250) 

0 

0 

111 

18816 

85 

0 

0 

D 

[3000,  5000) 

[-250,  250) 

0 

0 

1 

152 

1 

0 

0 

O 

[3000,  5000) 

[-250,  250) 

0 

1 

14992 

5695017 

8068 

0 

0 

B 

[5000,  oo] 

[-250,  250) 

0 

0 

33 

12535 

36 

0 

0 

C 

[5000,  oo] 

[-250,  250) 

0 

0 

0 

0 

0 

0 

0 

D 

[5000,  oo] 

[-250,  250) 

0 

0 

0 

0 

0 

0 

0 

O 

[5000,  oo] 

[-250,  250) 

0 

0 

5706 

2621051 

2624 

2 

0 

B 

[500,  1200) 

[250,  750) 

0 

0 

0 

479 

23748 

96 

0 

C 

[500,  1200) 

[250,  750) 

0 

0 

0 

1153 

70327 

147 

0 

D 

[500,  1200) 

[250,  750) 

0 

0 

0 

38448 

1993067 

8258 

0 

O 

[500,  1200) 

[250,  750) 

0 

0 

0 

59198 

3089682 

11222 

0 

B 

[1200,  3000) 

[250,  750) 

0 

0 

0 

202 

12523 

39 

0 

C 

[1200,  3000) 

[250,  750) 

0 

0 

0 

580 

39857 

177 

0 

D 

[1200,  3000) 

[250,  750) 

0 

0 

0 

6300 

345205 

1905 

0 

O 

[1200,  3000) 

[250,  750) 

0 

0 

0 

36232 

1845196 

8232 

2 

B 

[3000,  5000) 

[250,  750) 

0 

0 

0 

185 

8620 

63 

0 

C 

[3000,  5000) 

[250,  750) 

0 

0 

0 

96 

2607 

23 

0 

D 

[3000,  5000) 

[250,  750) 

0 

0 

0 

1 

37 

1 

0 

O 

[3000,  5000) 

[250,  750) 

0 

0 

0 

9269 

462258 

2288 

2 

B 

[5000,  oo] 

[250,  750) 

0 

0 

0 

42 

2273 

13 

0 

C 

[5000,  oo] 

[250,  750) 

0 

0 

0 

0 

0 

0 

0 

D 

[5000,  oo] 

[250,  750) 

0 

0 

0 

0 

0 

0 

0 

O 

[5000,  oo] 

[250,  750) 

0 

0 

0 

3037 

137408 

988 

2 

B 

[500,  1200) 

[750,  1250) 

0 

0 

0 

0 

162 

3778 

25 

C 

[500,  1200) 

[750,  1250) 

0 

0 

0 

0 

229 

7429 

31 

D 

[500,  1200) 

[750,  1250) 

0 

0 

0 

0 

10921 

304726 

1171 

O 

[500,  1200) 

[750,  1250) 

0 

0 

0 

1 

13859 

388637 

1555 

B 

[1200.  3000) 

[750,  1250) 

0 

0 

0 

0 

40 

1464 

12 

C 

[1200,  3000) 

[750,  1250) 

0 

0 

0 

0 

188 

6213 

38 

D 

[1200,  3000) 

[750,  1250) 

0 

0 

0 

0 

2059 

68013 

462 

O 

[1200,  3000) 

[750,  1250) 

0 

0 

0 

1 

8705 

256042 

1900 

B 

[3000,  5000) 

[750,  1250) 

0 

0 

0 

0 

68 

1059 

26 

C 

[3000,  5000) 

[750,  1250) 

0 

0 

0 

0 

25 

360 

10 

D 

[3000,  5000) 

[750,  1250) 

0 

0 

0 

0 

1 

9 

1 

Continued  on  next  page.  .  . 
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TABLE  B-3.  Continued 


A 

L 

h(t) 

[-2000,  -1250) 

[-1250,  -750) 

[-750,  -250)  [-250,  250) 

[250,  750) 

[750,  1250) 

[1250,  2000] 

O 

[3000,  5000) 

[750,  1250) 

0 

0 

0 

0 

2423 

68127 

717 

B 

[5000,  oo ] 

[750,  1250) 

0 

0 

0 

0 

16 

258 

5 

C 

[5000,  oo] 

[750,  1250) 

0 

0 

0 

0 

0 

0 

0 

D 

[5000,  oo] 

[750,  1250) 

0 

0 

0 

0 

0 

0 

0 

O 

[5000,  oo] 

[750,  1250) 

0 

0 

0 

1 

1105 

31246 

544 

B 

[500,  1200) 

[1250,  2000] 

0 

0 

0 

0 

0 

44 

2282 

C 

[500,  1200) 

[1250,  2000] 

0 

0 

0 

0 

0 

73 

2046 

D 

[500,  1200) 

[1250.  2000] 

0 

0 

0 

0 

3 

1694 

42489 

O 

[500,  1200) 

[1250,  2000] 

0 

0 

0 

0 

3 

2151 

57621 

B 

[1200,  3000) 

[1250,  2000] 

0 

0 

0 

0 

0 

12 

556 

C 

[1200.  3000) 

[1250.  2000] 

0 

0 

0 

0 

0 

36 

1177 

D 

[1200.  3000) 

[1250.  2000] 

0 

0 

0 

0 

0 

503 

15598 

O 

[1200,  3000) 

[1250,  2000] 

0 

0 

0 

0 

2 

1980 

63193 

B 

[3000,  5000) 

[1250.  2000] 

0 

0 

0 

0 

0 

29 

986 

C 

[3000,  5000) 

[1250,  2000] 

0 

0 

0 

0 

0 

10 

482 

D 

[3000,  5000) 

[1250,  2000] 

0 

0 

0 

0 

0 

1 

15 

O 

[3000.  5000) 

[1250,  2000] 

0 

0 

0 

0 

1 

733 

27729 

B 

[5000,  oo] 

[1250,  2000] 

0 

0 

0 

0 

0 

8 

252 

C 

[5000,  oo] 

[1250,  2000] 

0 

0 

0 

0 

0 

0 

0 

D 

[5000,  oo] 

[1250,  2000] 

0 

0 

0 

0 

0 

0 

0 

O 

[5000,  oo] 

[1250,  2000] 

0 

0 

0 

0 

0 

577 

24900 

We  look  at  the  row  of  interest  and  randomly  assign  a  value  to  h(t  +  1)  with  probability 
proportional  to  the  corresponding  elements  plus  1.  We  then  assign  xp(t  -hi)  to  a  random  value 
conditional  on  A,  xp{t),  L ,  and  h(t  -hi).  Finally,  we  assign  v(t  -h  1)  based  on  the  assignments  of  L, 
v(t),  h(t+  1),  and  ip(t- h  1).  The  process  is  repeated  for  each  time  step.  The  length  of  the  trajectory 
may  be  made  as  long  as  desired. 

Sampling  from  the  Bayesian  networks  produces  a  sequence  of  assignments  of  variables  to 
bins,  e.g.  [250,750).  As  described  above,  the  next  step  is  to  sample  uniformly  within  the  bins  to 
produce  real  values  (with  the  exception  of  a  bin  spanning  0  in  which  case  0  itself  is  selected).  If  we 
simply  sample  within  each  bin  at  each  time  step,  there  would  be  excessive  variability  in  the  vertical 
rates,  turn  rates,  and  acceleration.  We  therefore  only  resample  within  bins  at  rates  specified  in  the 
parameters  file: 

000  0.0211495  0.0412835  0.0885336 

The  first  three  variables  have  zero  as  their  rates  because  they  do  not  change  during  the  course  of 
the  trajectory.  The  last  three  elements  specify  the  rates  with  which  fe,  and  xp  change  within  bins. 
As  the  trajectory  evolves,  whenever  a  variable  switches  to  a  new  bin  we  must  sample  within  the 
bin.  Typically,  the  values  of  variables  continue  in  the  same  bin  at  the  next  time  step.  If  the  bin 
remains  the  same,  we  either  continue  with  the  same  value  or  resample  within  the  bin  according  to 
the  specified  rates. 

Once  the  initial  conditions  and  a  series  of  control  variables  (t>,  h,  and  ip)  have  been  selected, 
the  aircraft  trajectory  is  constructed  or  simulated  using  an  appropriate  dynamic  model.  The  control 
variables  are  assumed  to  be  held  constant  across  each  1  s  time  step. 
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APPENDIX  C 
BAYESIAN  NETWORKS 


This  appendix  briefly  reviews  Bayesian  networks.  Further  discussion  of  Bayesian  networks 
may  be  found  elsewhere  [19-21]. 

C.l  DEFINITION 

A  Bayesian  network  is  a  graphical  representation  of  a  multivariate  probability  distribution 
over  variables  X  =  X\, . . .  ,Xn.  In  particular,  a  Bayesian  network  is  a  directed  acyclic  graph  G 
whose  nodes  correspond  to  variables  and  edges  correspond  to  probabilistic  dependencies  between 
them.  Associated  with  each  variable  Xi  is  a  conditional  probability  distribution  P(xj  |  7T,),  where 
7 r7  denotes  an  instantiation  of  the  parents  of  Xi  in  the  graph.  The  probability  of  an  instantiation 
of  the  variables  is  specified  directly  by  the  conditional  probability  distributions  in  the  Bayesian 
network: 

n 

P(x)  =  P(x  1,...,X„)  =  I  7Ti).  (C-l) 

i=  1 


C.2  SAMPLING 

It  is  rather  straightforward  to  sample  from  a  multivariate  distribution  represented  by  a 
Bayesian  network.  The  first  step  is  to  produce  a  topological  sort  of  the  nodes  in  the  network. 
A  topological  sort  orders  the  nodes  in  a  Bayesian  network  such  that  if  a  node  Xj  comes  before 
Xj  there  does  not  exist  a  directed  path  from  Xj  to  Xi.  Every  Bayesian  network  has  at  least  one 
topological  sort,  but  there  may  be  many.  Efficient  algorithms  exist  for  finding  a  valid  topological 
sort  [22]. 

To  produce  a  sample  from  the  joint  distribution  represented  by  a  Bayesian  network,  we  simply 
iterate  through  a  topologically  sorted  sequence  of  the  variables  and  sample  from  their  conditional 
probability  distributions.  The  topological  sort  ensures  that  when  sampling  from  each  conditional 
probability  distribution  the  necessary  parents  have  been  instantiated. 


C.3  PARAMETER  LEARNING 

The  parameters  9  of  a  Bayesian  network  determine  the  associated  conditional  probability 
distributions.  Given  some  fixed  network  structure  G,  we  can  learn  these  parameters  from  data.  In 
this  appendix,  we  assume  that  the  variables  are  discrete. 

Before  discussing  how  to  learn  the  parameters  of  a  Bayesian  network,  it  is  necessary  to 
introduce  some  notation.  Let  r7  represent  the  number  of  instantiations  of  Xi  and  q7  represent 
the  number  of  instantiations  of  the  parents  of  Xi.  If  Xi  has  no  parents,  then  <?7  =  1.  The  j th 
instantiation  of  the  parents  of  Xi  is  denoted  7T 7j. 
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There  are  riQi  parameters  in  a  Bayesian  network.  Each  parameter  is  written  Ojjk  and 
determines  P(Xx  =  k  |  7r 7J),  i.e., 

P{Xi  =  k  |  7 Tij)  =  9ijk  . 

Although  there  are  parameters,  only  Y17=i(ri  ~  l)<ft  are  independent. 

Computing  the  posterior  p(6  \  D ,  G)  involves  specifying  a  prior  p(6  \  G )  and  applying  Bayes’ 


rule 


p(0\D,G)  = 


P(D  |  0,G)p(0  |  G) 
P{D  |  G) 


P(D  |  O.G)p{6  |  G) 
f  P(D  |  0,G)p(0  |  G)d0  ' 


(C-2) 


If  Nijk  is  the  count  of  Xi  =  A:  given  7T?j  in  the  data  D ,  then  the  probability  of  the  data  given 
the  parameters  6  is 


n  qi 


p(°  i 0) = nn  n  °Nijk 


ijk 


i=lj=lk=l 


(C-3) 


Let  Oij  =  (diji1 . . . ,  9{jri ) .  Since  6^  is  independent  of  when  ij  ^  i'j',  the  prior  probability 
of  the  parameters  assuming  a  fixed  structure  G  is 


n  qi 


p(0|G)  =  nn^;lG)- 

i=lj=l 


(C-4) 


The  density  p(0?j  |  G)  is  a  distribution  over  relative  frequencies.  Under  some  very  weak  assump¬ 
tions,  it  is  possible  to  prove  that  p{6Xj  \  G)  is  Dirichlet  (see  [21],  Section  6.2.3).  Hence, 


p(Oij  I  G)  = 


Ua,7o) 

rUTTn^uT) 


o 


if  0  <  dljk  <  1  and  Oijk  =  1 

otherwise 


where  a^i, . . . ,  aijTi  are  the  parameters  of  the  Dirichlet  distribution  and  axjo  =  Y7k=iaijk-  For 
the  prior  to  be  objective  (or  noninformative),  the  parameters  axjk  must  be  identical  for  all  k. 
Different  objective  priors  have  been  used  in  the  literature.  Cooper  and  Herskovits  [18]  use  =  1. 
Heckerman,  Geiger,  and  Chickering  [23]  use  and  justify  atjk  =  l/(r^j). 

It  is  possible  to  show  that  p(07J  |  D ,  G )  is  Dirichlet  with  parameters  &ijk+Nijk, . . . ,  a^+A^.. 
Hence, 


p(Oij\D,G)  = 


)Q  +  ATjj )  T  TT*i  nQijk^~^ijk  1 

ll*=i  r(«ofc+ATi>fc)  1  U=1  °ijk 
0 


if  0  <  0ijk  <  1  and  Y%=i  hjk  =  1 
otherwise 


where  Ntj  =  Ylk=i  Nijk. 

Sampling  from  a  Bayesian  network  with  G  known,  0  unknown,  and  D  observed  involves 
assigning  k  to  Xj  with  probability 


P{Xi  =  k  I  7 Tij,  D,G)  =  J  0ijkp(0ijk  I  D,  G)  dQijk  = 


°Hjk  +  Njjk 


^2k'=l(aijk'  +  Nijk') 


(C-5) 
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C.4  STRUCTURE  LEARNING 


Finding  the  most  likely  structure  G  that  generated  a  set  of  data  D.  The  objective  is  to  find 
the  most  likely  graph  given  data.  By  Bayes’  rule, 

P{G  |  D)  a  P{G)P(D  |  G)  =  P{G)  j  P(D  \  0.  G)p{6  \  G)  <10  .  (C-6) 

The  previous  section  explains  how  to  compute  the  likelihood  P(D  |  0,G)  and  the  prior  p(0  \  G). 
Cooper  and  Herskovits  [18]  show  how  to  evaluate  the  integral  above,  resulting  in 


P(G 


n  qt 


£)  =  P(G)Iin 

i= 1  j= 1 


r(c*ijo)  t~t  r (oLjjk  +  Nijk) 
r(a^o  +  Nij)  AA  r(aij^) 


(C-7) 


where  NtJ  =  Nijk-  Heckerman,  Geiger,  and  Chickering  [23]  suggest  priors  over  graphs,  but 

it  is  not  uncommon  in  the  literature  to  assume  a  uniform  prior.  For  numerical  convenience,  most 
Bayesian  network  learning  packages  calculate  and  report  log  P(G  \  D )  4-  A",  where  K  is  a  constant 
independent  of  G.  This  quantity  is  often  called  the  Bayesian  score  and  may  be  used  for  structure 
comparison  and  search. 


47 


APPENDIX  D 
DENSITY  PLOTS 


The  figures  in  this  appendix  show  the  density  of  aircraft  squawking  VFR  in  the  radar  data 
used  to  construct  our  model.  Density  is  measured  in  number  of  aircraft  per  NM3.  Shown  in  the 
plots  are  the  base-10  logarithm  of  the  density.  The  logarithm  of  the  density  was  used  because  of 
the  extreme  concentrations  in  some  areas  of  the  country,  particularly  Florida,  in  and  around  Dallas 
and  Houston,  Texas,  and  southern  California.  Each  pixel  measures  1/6  of  a  degree  on  a  side  in 
latitude  and  in  longitude.  Depending  on  latitude,  the  square  mileage  of  a  pixel  may  be  between 
64  NM2  at  the  northernmost  latitude  to  92  NM2  at  the  southernmost  latitude.  The  highest  density 
pixel  is  in  Miami-Dade  County  and  has  a  density  of  approximately  10“2  6  aircraft  per  NM3  over  all 
altitude  layers  between  500  ft  AGL  and  FL180,  with  the  highest  density  layer  between  500  ft  AGL 
and  1200  ft  AGL. 

It  is  important  to  emphasize  that  this  density  map  is  the  observed  density  from  our  data  set. 
Factors  such  as  terrain  and  curvature  of  the  earth  prevent  the  observation  of  aircraft,  particularly 
low  flying  aircraft  such  as  those  flying  VFR.  The  locations  of  the  sensors  are  indicated  on  the  maps 
by  magenta  circles.  It  is  no  coincidence  that  the  measured  density  in  the  lowest  altitude  layer,  for 
example,  tend  to  be  concentrated  around  the  locations  of  the  sensors.  In  particular,  coverage  at  low 
altitudes  tends  to  be  sparser  in  the  western  portion  of  the  country.  Coverage  above  5000  ft  AGL, 
however,  is  fairly  good  across  the  country. 

Certain  areas  on  the  maps  are  worth  attention.  Due  to  the  protected  area  around  Washington 
DC,  there  are  virtually  no  1200  squawking  aircraft  at  any  altitude  layer  in  that  region.  By  contrast, 
there  are  is  a  large  concentration  in  Florida  and  southern  Arizona,  especially  near  the  two  locations 
of  the  Einbry-Riddle  flight  school  in  Daytona  Beach,  Florida  and  Prescott,  Arizona.  In  addition, 
the  effect  of  the  “upside  down  wedding-cake”  airspace  class  structure  around  major  airports  can 
be  seen  if  one  compares  the  heavy  concentration  near  the  pixel  representing  Atlanta  and  Dallas 
at  lower  altitude  layers  against  the  empty  airspace  in  the  same  pixel  at  the  highest  altitude  layer. 
Finally,  it  can  be  seen  that  the  density  near  Denver  at  higher  altitudes  AGL  is  quite  low;  this  may 
be  attributed  to  the  fact  that  an  additional  5,000  ft  above  ground  in  that  region  puts  the  aircraft 
above  10,000  ft  MSL. 
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FIGURE  D-l.  VFR  density  over  all  altitude  layers  up  to  FL180  in  log10  ai 
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FIGURE  D-2.  VFR  density  between  500  and  1200ftAGL  in  log10  aircraft  per  NM'3. 
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FIGURE  D-3.  VFR  density  between  1200  and  3000ftAGL  in  log10  aircraft  per  NM3. 
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FIGURE  D-4.  VFR  density  between  3000  and  5000ft  AGL  in  log10  aircraft  per  NM3. 
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FIGURE  D-5.  VFR  density  between  5000ft  AG L  and  10,000ft  MSL  in  log10  aircraft  per  NM3. 
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FIGURE  D-6.  VFR  density  between  10,000 ft  MSL  and  FL180  in  log10  aircraft  per  NM3. 


APPENDIX  E 
SENSORS 


The  following  table  lists  the  RADES  sensors  that  contributed  to  our  model.  Shown  are  the 
number  of  hours  of  flight  time  by  transponder-equipped  aircraft  squawking  1200  captured  by  each 
sensor  between  December  1,  2007  and  December  7,  2007. 


TABLE  E-l.  Beacon  Hours 

60  NM  range  sensors 

628.19  □ 

82.64  D 

26.23 

426.58  C 

96.52  0 

291.51  □ 

295.06  □ 

165.55  0 

56.54  1 

72.57  0 

322.48  □ 

795.56  ■=□ 

461.82  □ 

977.69  BZD 

218.21  D 

236.84  D 

913.5  ■=□ 

1994.23  mz: . H 

247.43  D 

454.24  □ 

684.51  ED 

183.77  0 

201.24  0 

67.78  0 

67.23  ! 

206.52  0 

14.02 

224.52  D 

23.99 

201.91  □ 

ACT  Waco,  TX,  USA  (ASR-11) 

ADW  Camp  Springs  (Andrews  AFB),  MD,  USA  (ASR-9) 

BGR  Bangor,  ME,  USA  (ASR-8) 

BOS  Logan  Inti  (Boston),  MA,  USA  (ASR-9) 

BUF  Buffalo  Inti,  NY,  USA  (ASR-9) 

BWI  Baltimore  (BWI),  MD,  USA  (ASR-9) 

COS  Colorado  Springs,  CO,  USA  (ASR-11) 

COU  Columbia  MO,  MO,  USA  (ASR-11) 

CRW  Charleston,  WV,  USA  (ASR-8) 

DCA  Washington  National,  DC,  USA  (ASR-9) 

DOV  Dover  AFB,  DE,  USA  (GPN-20) 

EWR  Newark,  NJ,  USA  (ASR-9) 

GRK  Ft  Hood,  TX,  USA  (ASR-9) 

HPN  White  Plains  (Westchester  Co),  NY,  USA  (ASR-9) 

HUF  Terre  Haute,  IN,  USA  (ASR-8) 

IAD  Washington  Dulles  (Chantilly),  DC,  USA  (ASR-9) 

JFK  New  York  (JFK),  NY,  USA  (ASR-9) 

LAXB  Los  Angeles  Inti  -  North,  CA,  USA  (ASR-9) 

LFT  Lafayette,  LA,  USA  (ASR-11) 

MDT  Harrisburg,  PA,  USA  (ASR-9) 

MHT  Manchester,  NH,  USA  (ASR-9) 

MLU  Monroe  LA,  LA,  USA  (ASR-8) 

MRB  Martinsburg,  WV,  USA  (ASR-9) 

MUO  Mountain  Home  AFB,  ID,  USA  (GPN-20) 

NQXA  Key  West  NAS,  FL,  USA  (GPN-27) 

PWM  Portland  (Cumberland),  ME,  USA  (ASR-9) 

RCA  Ellsworth  AFB  (Rapid  City),  SD.  USA  (GPN-20) 

SAV  Savannah,  GA,  USA  (ASR-8) 

SAW  Marquette  (KI  Sawyer),  MI,  USA  (ASR-7F) 

SDF  Louisville,  KY,  USA  (ASR-9) 

Continued  on  next  page.  .  . 
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TABLE  E-l.  Continued 

229.89  0 

250.74  D 

414.31  □ 

SJU  San  Juan,  PR,  USA  (ASR-8) 

STT  St. Thomas,  PR,  USA  (ASR-8) 

SWF  Newburgh  (Stewart),  NY,  USA  (ASR-9) 

200  NM  range  sensors 

289.27  □ 

2274.9  K_,.  . 1 

92.82  0 

267.37  D 

133.24  D 

113.89  0 

703.96  EZI 

2584.5  Ml.  1 

403.56  O 

1335.22  gT~1 

479.19  O 

121.43  0 

744.49  EZI 

108.98  0 

744.27  EZI 

418.95  □ 

1605.12  1 

1282.42  CT1 

831.23  EZI 

503.78  C 

273.38  D 

199.97  D 

605.39  CD 

595.5  EZ 

893.41  KU 

201.11  D 

511.27  O 

82.66  0 

501.35  C 

154.15  D 

461.99  □ 

880.29  MZD 

468.43  O 

751.66  ■=□ 

769.26  EZ 

AEX  Alexandria,  LA,  USA  (AN/FPS-20A) 

ATL  Atlanta  (Marietta),  GA,  USA  (ARSR-1E) 

BAM  Battle  Mountain,  NV,  USA  (ARSR-2) 

CPV  Coopersville,  MI,  USA  (AN/FPS-66A) 

DSV  Buffalo  (Dansville),  NY,  USA  (ARSR-1E) 

FLX  Fallon,  NV,  USA  (AN/FPS-66A) 

FPK  Salt  Lake  City  (Francis  Peak),  UT,  USA  (ARSR-1E) 
FTW  Ft  Worth  (Keller),  TX,  USA  (ARSR-1E) 

GUP  Gallup  (Farmington),  NM,  USA  (ARSR-2) 

HOU  Houston  (Ellington),  TX,  USA  (ARSR-1E) 

IND  Indianapolis,  IN,  USA  (ARSR-1E) 

IRK  Kirksville,  MO,  USA  (ARSR-3) 

JOL  Elwood  (Joliet),  IL,  USA  (ARSR-3) 

LSK  Lusk,  WY,  USA  (ARSR-2) 

MGM  Montgomery,  AL,  USA  (ARSR-1E) 

PIT  Pittsburgh  (Oakdale),  PA,  USA  (AN/FPS-67B) 

QAS  Las  Vegas  (Angel  Peak),  NV,  USA  (AN/FPS-20E) 
QBE  Roanoke  (Bedford),  VA,  USA  (ARSR-3) 

QBN  Binns  Hall,  VA,  USA  (ARSR-3) 

QBZ  Oskaloosa,  KS,  USA  (ARSR-2) 

QCF  Du  Bois  (Clearfield),  PA,  USA  (ARSR-3) 

QCK  Boise  (Cascade),  ID,  USA  (ARSR-2) 

QDB  Cleveland  (Brecksville),  OH,  USA  (ARSR-1E) 

QDT  Detroit  (Canton),  MI,  USA  (ARSR-1E) 

QHA  Cummington,  MA,  USA  (AN/FPS-67B) 

QHB  St.  Albans,  VT,  USA  (AN/FPS-67B) 

QHZ  Horicon,  WI,  USA  (ARSR-2) 

QJB  Gettysburg,  SD,  USA  (AN/FPS-67B) 

QJE  Minneapolis  (Apple  Valley),  MN,  USA  (ARSR-1E) 

QJO  Arlington,  IA,  USA  (ARSR-3) 

QJQ  San  Juan  (Pico  del  Este),  PR,  USA  (AN/FPS-20E) 
QNK  Lincolnton,  GA,  USA  (ARSR-3) 

QNM  Newport,  MS,  USA  (ARSR-3) 

QOJ  Nashville  (Joelton),  TN,  USA  (ARSR-1E) 

QPC  Haleyville,  AL,  USA  (AN/FPS-67B) 

Continued  on  next  page.  . . 

58 


1119.91  M — I 
903.24  MU 
1547.86  1.  1 

346.34  □ 
681.5  MU 
823.63  MU 
548.93  IZ 
388.89  □ 
97.39  0 
510.95  IZ 
221.78  D 
146.17  0 
613.66  IZ) 
430.32  E 
334.63  □ 
373.29  □ 
1401.67  1 

767.89  EZ 
735.82  MU 
265.31  □ 


1320.68 


5.27 


1220.64  MUD 
350.53  □ 
317.16  D 
3833.04  T"  . . "'1 

566.36  IZ 
432.44  O 
2632.79  Mk  I 

662.4  MU 
942.28  MU 
689.27  MU 
4235.19  I 


1321.28  g - 1 

112.12  D 
125.28  D 
159.75  0 


TABLE  E-l.  Continued 

QPK  Denver  (Parker),  CO,  USA  (ARSR-1E) 

QPL  The  Plains,  VA,  USA  (ARSR-3) 

QRB  Citronelle,  AL,  USA  (ARSR-2) 

QRC  Benton,  PA,  USA  (AN/FPS-67B) 

QRI  Lynch,  KY,  USA  (ARSR-2) 

QRL  Raleigh  (Benson),  NC,  USA  (ARSR-1E) 

QRM  Charlotte  (Maiden),  NC,  USA  (ARSR-1E) 

QSA  Albuquerque  (West  Mesa),  NM,  USA  (AN/FPS-66A) 
QSI  Lovell,  WY,  USA  (ARSR-2) 

QTZ  La  Grange,  IN,  USA  (ARSR-1E) 

QUZ  Hanna  City,  IL,  USA  (AN/FPS-67B) 

QWC  Mesa  Rica,  NM,  USA  (ARSR-1E) 

QWO  London,  OH,  USA  (ARSR-1E) 

QXR  Russellville,  AR,  USA  (AN/FPS-64A) 

QXS  Odessa  (Andrews),  TX,  USA  (ARSR-1E) 

QYB  Byhalia  (Memphis),  MS,  USA  (ARSR-1E) 

QYS  Rogers,  TX,  USA  (ARSR-1E) 

SEA  Seattle  (Ft  Lawton),  WA,  USA  (ARSR-1E) 

STL  St,  Louis,  MO,  USA  (ARSR-1E) 

SVC  Silver  City,  NM,  USA  (ARSR-2) _ 

250  NM  range  sensors 
A  JO  Ajo,  AZ,  USA  (ARSR-4) 

BAR  Barrington,  NS,  Canada  (AN/FPS-117) 

CTY  Cross  City,  FL,  USA  (ARSR-4) 

DMN  Deming  (Magdelina  Peak),  NM,  USA  (ARSR-4) 

LCH  Lake  Charles,  LA,  USA  (ARSR-4) 

MLB  Melbourne,  FL,  USA  (ARSR-4) 

NEN  Jacksonville  (White  House  FI),  FL,  USA  (ARSR-4) 
NEW  New  Orleans  (Slidell),  LA,  USA  (ARSR-4) 

NQX  Key  West,  FL,  USA  (ARSR-4) 

NSD  San  Clemente  Island,  CA,  USA  (ARSR-4) 

PAM  Panama  City  (Tyndall),  FL,  USA  (ARSR-4) 

PRB  Paso  Robles,  CA,  USA  (ARSR-4) 

QEA  North  Truro,  MA,  USA  (ARSR-4) 

QFN  Ft  Green,  FL,  USA  (ARSR-4) 

QIE  Gibbsboro,  NJ,  USA  (ARSR-4) 

QJA  Empire,  MI,  USA  (ARSR-4) 

QJD  Nashwauk,  MN,  USA  (ARSR-4) 

QKW  Makah,  WA,  USA  (ARSR-4) _ 

Continued  on  next  page.  .  . 
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TABLE  E-l.  Continued 

3158.42  K _ 1 

2649.29  Ml  1 

692.15  EZ] 

221.85  D 

457.99  □ 

328.46  □ 

1129.22  H — 1 

1173.37  Ml — 1 

541.18  □ 

41.6 

54.44  1 

34.93 

20.77 

225.25  □ 

66.18  D 

198.62  □ 

337.83  □ 

QM8  Tamiami,  FL,  USA  (ARSR-4) 

QMV  Mill  Valley,  CA,  USA  (ARSR-4) 

QNA  Morales,  TX,  USA  (ARSR-4) 

QNW  Eagle  Peak,  TX,  USA  (ARSR-4) 

QOM  King  Mountain,  TX,  USA  (ARSR-4) 

QRJ  Charleston  (Jedburg),  SC,  USA  (ARSR-4) 

QRW  Mt.  Laguna,  CA,  USA  (ARSR-4) 

QVH  Riverhead,  NY,  USA  (ARSR-4) 

QVR  Oceana,  VA,  USA  (ARSR-4) 

QWA  Watford  City,  ND,  USA  (ARSR-4) 

QXU  Utica  (Remsen),  NY,  USA  (ARSR-4) 

QYA  Bucks  Harbor,  ME,  USA  (ARSR-4) 

QYD  Caribou,  ME,  USA  (ARSR-4) 

QZA  Oilton,  TX,  USA  (ARSR-4) 

QZZ  Rainbow  Ridge,  CA,  USA  (ARSR-4) 

RSG  Rocksprings,  TX,  USA  (ARSR-4) 

SLE  Salem  (Laurel  Mountain),  OR,  USA  (ARSR-4) 
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APPENDIX  F 
TRACKING  AND  FUSION 


Converting  raw  radar  reports  into  tracks  that  are  usable  for  our  model  development  is  a 
two-stage  process.  The  first  stage  involves  forming  local  tracks  from  the  reports  associated  with 
each  sensor.  The  second  stage  involves  fusing  the  local  tracks  from  multiple  sensors  to  form  global 
tracks.  This  appendix  provides  a  brief  overview  of  tracking  and  fusion. 

We  use  the  tracking  algorithms  from  the  Mode  S  and  ASR-9  systems  [13],  the  two  most 
modern  sensors  in  the  Air  Traffic  Control  System.  The  beacon  correlation  algorithms  come  from 
Mode  S  and  the  primary  radar  correlation  algorithms  come  from  ASR-9  Processor  Augmentation 
•  Card  (9-PAC).  Both  systems  are  integrated  to  provide  a  consistent  track,  although  for  the  purpose 

of  this  report  we  ignore  primary  only  reports. 

After  the  reports  of  each  sensor  have  been  correlated  into  local  tracks,  we  can  fuse  the  local 
tracks  to  provide  a  global  picture  of  the  airspace.  Fusion  performs  the  following  functions: 

1.  Merge  tracks  from  multiple  sensors  that  correspond  to  the  same  aircraft  into  a  single  global 
track. 

2.  Compute  the  speed  and  heading  of  each  global  track  to  permit  trajectory  predic  tions. 

3.  Correct  sensor  tracking  errors  that  would  have  led  to  split  global  tracks  and  thus  false  en¬ 
counters. 

We  use  a  track-to-track  fusion  method,  meaning  that  we  track  each  sensor’s  individual  reports 
and  then  we  merge  all  of  the  tracks  [24].  The  main  advantages  of  track-to-track  fusion  over  report- 
to-track  fusion,  in  which  all  reports  are  correlated  directly  to  global  tracks,  are: 

1.  Bias  independence  for  velocity  determination  and  maneuver  detection. 

2.  Removal  of  short  update  interval  velocity  anomalies. 

3.  Reduced  likelihood  of  forming  clutter  tracks. 

4.  Reduced  likelihood  of  introducing  incorrect  data  points  into  the  track  due  to  correlation 
errors. 


Track  merging  employs  position,  velocity,  Mode  A  code,  and  altitude  as  matching  attributes  over 
the  entire  track.  Track  merging  for  tracks  with  discrete  codes,  which  are  unique  within  an  area, 
employs  large  correlation  boxes  for  each  of  the  matching  attributes.  Other  tracks  (1200  code 
or  radar-only)  must  pass  more  stringent  position  tests  and  velocity  tests  in  order  to  be  merged 
together. 

Our  fusion  method  works  forward  in  time.  Thus,  there  is  often  doubt  in  whether  or  not 
tracks  from  multiple  radars  are  indeed  the  same  track  with  just  a  few  data  points.  In  cases  of 
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doubt,  tentative  matches  are  remembered  that  can  be  upgraded  to  a  merge  after  more  scans  of 
information  are  obtained.  Merges  are  checked  each  scan  and  can  be  undone  if  later  found  to  be 
unsatisfactory.  Aircraft  code  changes  are  also  accommodated,  although  they  must  be  verified  by 
other  sensors  in  the  merge  set  of  the  global  track  before  being  accepted.  If  only  one  sensor  reports 
a  code  change,  it  is  assumed  that  the  sensor  had  a  track  swap,  and  the  merge  situation  is  altered 
accordingly. 

The  remainder  of  this  section  explains  how  we  add  local  sensor  tracks  to  global  tracks,  how 
to  estimate  track  velocity  as  part  of  the  fusion  process,  and  how  to  filter  for  encounters. 


F.l  ADDING  LOCAL  SENSOR  TRACKS  TO  GLOBAL  TRACKS 

In  order  to  facilitate  fusing  of  tracks  from  multiple  sensors  into  global  tracks,  we  break  the 
continental  United  States  into  20  NM  by  20  NM  bins.  Every  track  is  associated  to  a  geographic 
bin.  Whenever  the  fusion  process  receives  a  new  local  sensor  track,  or  a  later  report  for  an  as  yet 
unfused  local  track,  an  attempt  is  made  to  fuse  this  track  to  an  existing  global  track.  This  process 
is  performed  by  comparing  the  new  track  to  all  neighboring  global  tracks.  The  neighboring  global 
tracks  for  a  local  track  include  all  global  tracks  in  the  same  geographical  bin  as  the  local  track 
plus  global  tracks  in  the  surrounding  bins  (totalling  9  bins).  Several  tests  must  be  passed  for  the 
successful  fusion  of  the  single  track  to  a  global  track: 

1.  The  global  track  must  not  already  be  fused  to  another  local  track  from  the  new  track’s  sensor. 

2.  The  tracks  must  agree  on  Mode  A  code  (primary-only  tracks  automatically  pass). 

3.  If  the  code  agreement  was  on  a  discrete  code,  a  very  coarse  horizontal  positional  test,  must 
be  passed. 

4.  If  the  code  agreement  was  on  1200  code,  or  no  codes,  a  tighter  positional  test  must  be  passed, 
as  well  as  altitude  and  velocity  tests.  However,  if  only  the  velocity  test  fails,  a  potential  fusion 
is  declared;  3  successive  potential  fusions  with  the  same  global  track  results  in  a  successful 
fusion. 

If  more  than  one  possible  global  track  satisfies  the  fusion  tests,  the  one  with  the  highest  matching 
score  is  chosen.  Existing  fusion  matches  are  checked  each  time  a  new  report  is  received.  If  the 
tests  fail  for  three  scans  in  a  row,  the  fusion  is  ended,  and  a  new  global  track  is  sought  for  the  local 
sensor  track. 


F.1.1  Code  Matching  Test 

Normally,  the  code  of  the  new  track  is  the  same  as  the  code  of  the  global  track.  However,  code 
mismatching  can  complicate  the  fusion  process.  Code  declaration  errors  due  to  data  corruption, 
missing  codes,  and  code  changes  due  to  controller  action  are  all  common.  For  this  reason,  associated 
with  each  global  track  is  an  established  code  and  an  alternate  code,  which  is  the  code  of  the  most 
recent  report.  Usually  these  two  codes  are  the  same.  When  an  alternate  code  is  different  from 
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the  established  code  for  three  successive  reports,  however,  the  established  code  is  updated  to  the 
alternative  code.  Reports  with  no  beacon  code  are  ignored  in  this  process.  However,  if  a  track  has 
never  had  been  associated  with  a  beacon  code  report,  we  consider  this  track  to  be  a  primary-only 
and  give  the  track  an  established  code  of  0. 

If  both  the  local  track  the  global  track  have  a  beacon  code,  then  a  successful  code  match  is 
declared  if  any  of  the  following  statements  are  true: 

.  1.  The  established  codes  match. 

2.  The  alternate  codes  match. 

#  3.  One  of  the  established  codes  and  the  other  alternate  code  match. 

However,  failure  is  declared  if  the  match  is  on  code  0,  and  both  tracks  have  a  beacon  code  in  the 
other  code  slot  that  do  not  match.  Lastly,  if  one  track  is  radar  only,  and  the  other  track  has  a 
beacon  code,  failure  is  declared. 

To  handle  local  track  code  changes  due  to  a  track  swap  in  the  single  sensor  tracker,  confir¬ 
mation  of  the  code  change  by  the  global  track  is  required.  If  at  the  time  of  the  local  track  code 
change,  the  global  track  has  had  an  update  by  a  different  sensor’s  track  with  the  old  beacon  code, 
a  track  swap  is  declared,  and  the  local  track  is  removed  from  fusion  with  the  global  track.  The 
local  track  then  undergoes  a  new  fusion  process. 

F.1.2  Horizontal  Position  Matching  Test 

Horizontal  positional  matching  requires  agreement  between  the  global  track's  most  recent 
horizontal  position  xg,yg  and  the  new  local  track’s  horizontal  position  x\,y\  projected  back  to  the 
time  of  the  global  track.  This  test  is  simple  if  both  track’s  reports  contained  altitude.  If  a  radar 
report  contains  altitude  /?,  range  p,  azimuth  6 ,  then  the  track’s  horizontal  position  x,y  can  be 
empirically  determined. 

If,  however,  the  altitude  of  a  track  is  unknown,  then  the  tracks'  altitude  has  to  be  assumed 
(or  guessed)  in  order  to  derive  the  track’s  horizontal  position.  Simply  guessing  an  altitude  may 
produce  a  erroneous  x,y  position.  In  order  to  use  a  reasonable  altitude  value,  we  employ  the 
following  algorithm: 

1.  If  only  one  track  has  known  altitude,  then  we  convert  the  other  track’s  stored  p,  0  position 
to  x,  y  using  the  first  track’s  altitude. 

^  2.  If  neither  track  has  known  altitude  (which  is  always  true  for  a  primary-only  match),  we 

consider  all  altitudes  from  ONM  to  7NM  at  1NM  steps.  We  then  use  the  altitude  that 
produces  the  closest  positional  match  between  the  two  tracks.  While  a  smaller  step  size  may 
produce  more  accurate  estimates,  we  found  that  1 NM  is  sufficient  for  fusing  two  tracks. 

Horizontal  positional  agreement  is  declared  if  the  horizontal  distance  between  the  two  tracks 
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is  less  than  an  acceptable  value: 


\j(H~  Zl)2  +  (yg  -  yi)2  <  Armax  +  Sa^pg  +  3(Tazp/ ,  (F-l) 

where  we  use  Armax  =  20  NM  for  a  discrete  code  match  and  A rmax  =  1  NM  for  a  1200  code  or 
radar  match.  The  standard  deviation  for  horizontal  position  error  terms  <raz  account  for  positional 
errors  due  to  azimuth  noise  in  radar  measurements  (which  is  the  dominant  source  of  horizontal 
position  error).  We  model  the  standard  deviation  of  the  azimuth  noise  as  3  milliradians  for  the 
data  format  of  our  radar  feed. 

F.1.3  Altitude  Matching  Test 

Altitude  matching  requires  agreement  between  the  local  and  global  track  altitudes  when 
both  are  known.  Two  comparisons  are  tested;  the  success  of  either  test  results  in  a  match.  The 
comparison  test  is: 


SO 

< 

= 

*g-*i 

A  hi 

= 

\h%  - 

A/i, 

< 

A/imax 

A/i, 
A  t\ 

< 

A/lmax 

where  we 

use 

A/?max  =  ( 

Ft  and  A/?max  =  100  ft/s.  Since  the  altitude  of  a  track  can  significantly 
change  between  sequential  reports,  we  test  both  the  most  recent  and  previous  local  track  altitudes 
with  the  most  recent  global  track  altitude  update.  Only  one  of  the  local  track  altitudes  is  required 
to  pass  the  test. 


F.1.4  Velocity  Matching  Test 

Velocity  matching  requires  agreement  between  the  two  track  headings  jp  and  speeds  s  accord¬ 
ing  to  the  following  tests: 


iV’g  -  V'll 

< 

A^max 

|sg  -  Sll 

< 

Asmax 

-  < 

<  2 

2  “ 

si 

4 

where  we  use  At/;max  =  45°  and  Asmax  =  100  kt.  The  last  test  is  needed  for  slow  aircraft  and 
clutter  tracks,  to  prevent,  for  example,  speeds  of  20  and  110  kt  from  agreeing. 


F.2  DETERMINING  TRACK  AIRSPEED  AND  HEADING 

Determining  a  global  track’s  airspeed  and  heading  is  a  two  step  process.  First,  the  individual 
sensor  tracks  are  smoothed.  Second,  the  individual  track  are  averaged  using  relative  weights  that 
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account  for  sensor  update  times  and  the  quality  of  each  sensor’s  measurement.  We  apply  both 
alpha  smoothing  and  curve  fitting  to  determine  airspeed  and  heading,  depending  upon  the  track 
situation.  We  have  implemented  various  maneuver  detection  algorithms,  and  tracking  is  dependent 
upon  the  current  turn  rate  and  acceleration  states  of  the  track. 

F.2.1  Local  Track  Smoothing 

First,  we  require  that  the  track  has  moved  a  minimum  distance  for  it  to  be  considered.  If 
the  track  never  moves  more  than  1  NM,  then  the  track  is  thrown  out.  After  the  movement  test  is 
satisfied,  the  track’s  airspeed  and  heading  are  calculated  from  the  new  and  previous  positions.  We 
then  update  the  local  track’s  airspeed  and  heading  estimates  using  alpha  smoothing. 

First,  the  current  heading  estimate  xp^\  and  its  difference  from  the  previous  estimate  xp^~l\ 
are  given  by 

xj)^  =  atan2  (y^  — 

A =  x/j^  —  xp^~^ 


Next,  we  determine  the  current  turn  rate  state  Sy  of  the  track: 

2  if  Axp^  >  ^heading 

1  if  Axp^  >  Axpm[n 

-2  if  Axp^  <  -^heading  (F-2) 

-1  if  A x/jW  <  -A^min 
0  otherwise 

where  we  use  A0mjn  =  3°  and  ^heading  is  the  standard  deviation  of  the  heading  noise,  which  is 
calculated  from  the  standard  deviations  for  range  and  azimuth  noise  of  the  sensors.  Note  that  a 
positive  A0  value  corresponds  to  a  right  turn,  while  a  negative  value  corresponds  to  a  left  turn. 
We  then  use  to  determine  the  smoothing  value  alpha  a  in  Table  F-l  of  the  individual  tracks 
that  will  be  used  to  calculate  the  heading  of  the  global  track  at  the  current  time. 

The  new  track  heading  is  finally  given  by: 

xp(j)  =  'ijjti-i)  -j_  &  x  A'ip^ 

and  we  iterate  through  this  process  for  the  entire  track. 

The  process  to  estimate  airspeed  s  is  similar,  with  one  important  difference.  If  successive 
positions  are  simply  connected,  then  the  airspeed  estimates  will  always  be  too  high,  since  the 
aircraft  will  appear  to  uzig-zag”  along  the  track.  Thus,  only  the  projection  of  the  velocity  vector 
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TABLE  F-l.  Smoothing  values  depending  on  the  current  and  previous  turn  states. 


Previous  State 

Current  Turn  State 

Large 

Left 

Turn  (-2) 

Small 

Left 

Turn  (-1) 

No  Turn 
(0) 

Small 

Right 

Turn  (+1) 

Large 

Right 

Turn  (+2) 

Large  Left  (-2) 

0.7 

0.7 

0.4 

0.5 

0.5 

Small  Left  (-1) 

0.7 

0.4 

0.4 

0.5 

0.5 

No  Turn  (0) 

0.4 

0.4 

0.3 

0.4 

0.4 

Small  Right  (+1) 

0.5 

0.5 

0.4 

0.4 

0.7 

Large  Right  (+2) 

0.5 

0.5 

0.4 

0.7 

0.7 

onto  the  track’s  heading  vector  is  used  to  determine  the  track’s  airspeed: 

Jj)  -  cos  (  A^U)  \  x  J(xU)  -  xU-1))2  +  (yU)  -y(j-»)2 
s  -  cos  ^  2  J  x  y  1) 

A  —  s^~1^ 


We  then  determine  the  current  airspeed  acceleration  state  Ss  of  the  track  using  a  similar 
technique  as  we  did  for  turn  rate. 


2 

1 

-2 

-1 

0 


if  Asti)  >  crspeed 
if  Asti)  >  Asmin 

if  As^^  <  —  ^heading 
if  Asti)  <  — Asmin 
otherwise 


where  we  use  Asmin  =  18  kt  and  crspeed  is  the  standard  deviation  of  airspeed  error  due  to  noise 
in  range  and  azimuth  measurements  from  the  radar  sensors.  The  speed  smoothing  and  the  speed 
alpha  table  rules  are  the  same  as  for  the  heading  case. 


F.2.2  Global  Track  Smoothing 

In  order  to  determine  a  global  track’s  airspeed  and  heading  at  each  measurement,  we  use  a 
weighted  least  squares  estimation  approach.  In  this  section  we  describe  in  detail  the  approach  for 
determining  the  tracks  heading. 

First,  each  sensor’s  heading  estimate  is  assigned  a  weight  at  the  current  time  t ^  as  follows: 


w 


(c) 


r-l 

heading 


^max 


JJ)  ,  ,0-1) 

Li  ^ 


t(c) 
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where  fading  is  the  standard  deviation  of  the  heading  noise  and  £max  =  18  s  is  a  discounting 
factor  that  takes  into  account  the  time  difference  between  the  measurement  from  the  sensor  being 
considered  and  the  time  for  when  we  are  determining  heading.  The  time  corresponds  to  the 
time  of  the  next  closest  measurement  for  sensor  i  with  respect  to  the  current  time  that  we  are 
trying  to  determine  the  track’s  heading. 

Next,  we  determine  the  total  turn  state  score  for  the  track 


where  TV  is  the  number  of  sensors  supporting  the  track  and  Si  is  the  current  turn  state  value  for 
the  ith  sensor  defined  in  Equation  F-2.  If  the  turn  rate  score  is  less  than  TV,  then  we  consider  the 
track  to  be  non-turning  and  the  current  global  heading  is  simply  the  weighted  average  of  the  TV 
sensor  heading  estimates: 


(c) 

global 


x  w- 


En 
i=\w 


( C ) 


Otherwise,  if  the  turn  rate  score  is  greater  than  or  equal  to  TV,  then  we  consider  the  track 
to  be  in  a  turn.  In  this  case,  we  utilize  weighted  least  squares  estimated  to  determine  a  first-order 
relationship  between  time  and  heading. 

The  global  track  speed  calculation  is  identical  in  form  to  the  global  track  heading  calculation. 
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APPENDIX  G 
NETWORK  CANDIDATES 


This  section  displays  a  selection  of  various  network  structures  along  with  their  log-Bayesian 
score,  number  of  edges,  and  number  of  parameters.  Bayesian  model  selection  balances  the  com¬ 
plexity  of  the  model  according  to  the  amount  of  data  available. 


G.l  INITIAL  NETWORK  CANDIDATES 


Score  =  -500006790  (best),  Edges  =  15,  Parameters  =  31359 
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Score  =  -502636693,  Edges  =  14,  Parameters  =  9855 


Score  =  -502638225,  Edges  =  14,  Parameters  =  9855 
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Score  =  -503170515,  Edges  =  13,  Parameters  =  7167 


Score  =  -503711030,  Edges  =  12,  Parameters  =  7083 
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Score  =  -507340429,  Edges  =  10,  Parameters  =  1775 


Score  =  -507958719,  Edges  =  10,  Parameters  =  1791 
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>4 

L 

V 

V 

h 

Score  =  -531540309,  Edges  =  0,  Parameters  =  29 


G.2  TRANSITION  NETWORK  CANDIDATES 


Score  =  -13243212  (best),  Edges  =11,  Parameters  =  9296 


* 


Score  =  -13266356,  Edges  =  10,  Parameters  =  5936 
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Score  =  -13276011,  Edges  =  12,  Parameters  =  18592 


Score  =  -13277700,  Edges  =  10,  Parameters  =  11312 
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Score  =  -13340477,  Edges  =  12,  Parameters  =  42784 


Score  =  -13380276,  Edges  =  14,  Parameters  =  74368 
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Score  =  -13942906,  Edges  =  21,  Parameters  =  7651840 


j>{t+ 1) 


h(t  + 1 ) 


v(t  + 1) 


Score  =  -180576482,  Edges  =  0,  Parameters  =  16 
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